module sealw99_mod
!
! fil
!
!
! ds
!
!
! This code provides the core functionality of FMS longwave
! radiation. It is based on exchange method with prescribed
! coefficients embedded in the code.
!
!
!
!
! 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, NOTE, close_file
use time_manager_mod, only: time_type, operator(>=), get_time, &
operator(-), get_date
use constants_mod, only: constants_init, diffac, radcon_mks, &
SECONDS_PER_DAY, radcon
! radiation package shared modules:
use rad_utilities_mod, only: rad_utilities_init, Lw_control, &
cldrad_properties_type, &
cld_specification_type, &
lw_output_type, longwave_tables1_type, &
longwave_tables2_type, &
longwave_tables3_type, atmos_input_type,&
Cldrad_control, radiative_gases_type, &
aerosol_type, aerosol_properties_type, &
aerosol_diagnostics_type, &
optical_path_type, gas_tf_type, &
lw_table_type, Lw_parameters, &
lw_diagnostics_type, lw_clouds_type, &
locate_in_table, looktab, mass_1, &
temp_1, Rad_control
use longwave_params_mod, only: longwave_params_init, NBCO215, &
NBLY_CKD, NBLY_RSB, longwave_params_end
! radiation package modules:
use longwave_clouds_mod, only: longwave_clouds_init, cldtau, cloud, &
lw_clouds_dealloc, longwave_clouds_end, &
thickcld
use longwave_fluxes_mod, only: longwave_fluxes_ks, &
longwave_fluxes_init, &
longwave_fluxes_end, &
longwave_fluxes_k_down, &
longwave_fluxes_KE_KEp1, &
longwave_fluxes_diag, &
longwave_fluxes_sum
use longwave_tables_mod, only: longwave_tables_init, &
longwave_tables_end
use optical_path_mod, only: optical_path_setup, &
optical_path_init, &
optical_trans_funct_from_KS, &
optical_trans_funct_k_down, &
optical_trans_funct_KE, &
optical_trans_funct_diag, &
get_totvo2, get_totch2o, get_totch2obd,&
optical_dealloc, optical_path_end
use gas_tf_mod, only: co2coef, transcolrow, transcol, &
get_control_gas_tf, gas_tf_init, &
gas_tf_dealloc, gas_tf_end, &
trans_sfc, trans_nearby
use lw_gases_stdtf_mod, only: lw_gases_stdtf_init, cfc_indx8, &
cfc_indx8_part, cfc_exact, &
lw_gases_stdtf_time_vary, &
lw_gases_stdtf_dealloc, co2_lblinterp, &
ch4_lblinterp, n2o_lblinterp, &
lw_gases_stdtf_end
!------------------------------------------------------------------
implicit none
private
!-------------------------------------------------------------------
! sealw99_mod is the internal driver for the 1999 sea longwave
! radiation code.
!---------------------------------------------------------------------
!---------------------------------------------------------------------
!----------- version number for this module -------------------
character(len=128) :: version = '$Id: sealw99.F90,v 19.0 2012/01/06 20:23:35 fms Exp $'
character(len=128) :: tagname = '$Name: tikal $'
logical :: module_is_initialized = .false.
!---------------------------------------------------------------------
!------- interfaces --------
public &
sealw99_init, sealw99_time_vary, sealw99, &
sealw99_endts, sealw99_end
private &
check_tf_interval, obtain_gas_tfs, &
sealw99_alloc, &
cool_to_space_approx, cool_to_space_exact, &
e1e290, e290, esfc, enear, &
co2_source_calc, nlte, co2curt, &
co2_time_vary, ch4_time_vary, n2o_time_vary
!---------------------------------------------------------------------
!-------- namelist ---------
logical :: &
do_thick = .false. ! perform "pseudo-convective
! adjustment" for maximally
! overlapped clouds ?
logical :: &
do_lwcldemiss = .false. ! use multiple bands to calculate
! lw cloud emissivites ?
logical :: &
do_ch4lbltmpint = .false. ! perform and save intermediate
! flux calculations for ch4?
logical :: &
do_n2olbltmpint = .false. ! perform and save intermediate
! flux calculations for n2o?
logical :: &
do_nlte = .false. ! there is a non-local thermodynamic
! equilibrium region at top
! of model ?
character(len=16) :: &
continuum_form = ' ' ! continuum specification; either
! 'ckd2.1', 'ckd2.4', 'mt_ckd1.0',
! 'rsb' or 'none'
character(len=16) :: &
linecatalog_form = ' ' ! line catalog specification; either
! 'hitran_1992' or 'hitran_2000'
real :: &
co2_tf_calc_intrvl = 1.0E6 ! interval between recalculating co2
! transmission functions, relevant
! for time-varying co2 cases
! [ hours ]
logical :: &
calc_co2_tfs_on_first_step = .true.
! always calculate co2 tfs on
! first time step of job ?
logical :: &
use_current_co2_for_tf = .false.
! use current co2 mixing ratio for
! calculating tfs ?
real :: &
co2_tf_time_displacement = 0.0
! time displacement from job start
! to the point in time where co2
! tfs are to be valid -- may be (-),
! 0.0 or (+); used only when
! calc_co2_tfs_on_first_step is true
! [ hours ]
real :: &
ch4_tf_calc_intrvl = 1.0E6 ! interval between recalculating ch4
! transmission functions, relevant
! for time-varying ch4 cases
! [ hours ]
logical :: &
calc_ch4_tfs_on_first_step = .true.
! always calculate ch4 tfs on
! first time step of job ?
logical :: &
use_current_ch4_for_tf = .false.
! use current ch4 mixing ratio for
! calculating tfs ?
real :: &
ch4_tf_time_displacement = 0.0
! time displacement from job start
! to the point in time where ch4
! tfs are to be valid -- may be (-),
! 0.0 or (+); used only when
! calc_ch4_tfs_on_first_step is true
! [ hours ]
real :: &
n2o_tf_calc_intrvl = 1.0E6 ! interval between recalculating n2o
! transmission functions, relevant
! for time-varying n2o cases
! [ hours ]
logical :: &
calc_n2o_tfs_on_first_step = .true.
! always calculate n2o tfs on
! first time step of job ?
logical :: &
use_current_n2o_for_tf = .false.
! use current n2o mixing ratio for
! calculating tfs ?
real :: &
n2o_tf_time_displacement = 0.0
! time displacement from job start
! to the point in time where n2o
! tfs are to be valid -- may be (-),
! 0.0 or (+); used only when
! calc_n2o_tfs_on_first_step is true
! [ hours ]
integer :: &
verbose = 0 ! verbosity level, ranges from 0
! (min output) to 5 (max output)
logical :: &
calc_co2_tfs_monthly = .false.
logical :: &
calc_ch4_tfs_monthly = .false.
logical :: &
calc_n2o_tfs_monthly = .false.
integer :: &
no_h2o_bands_1200_1400 = 1 ! number of bands in the lw par-
! ameterization between 1200 and 1400
! cm (-1); 0 and 1 have been
! tested. other potentially available
! values are 2, 4, 10 and 20.
logical :: &
use_bnd1_cldtf_for_h2o_bands = .false. ! the 1200-1400 cm(-1) band
! uses the same radiative
! properties as the 0-160,
! 1400-2200 band. needed
! for backward compatibil-
! ity
namelist / sealw99_nml / &
do_thick, do_lwcldemiss, &
! do_nlte, do_ch4n2olbltmpint, &
do_nlte, do_ch4lbltmpint, do_n2olbltmpint, &
continuum_form, linecatalog_form, &
verbose, &
no_h2o_bands_1200_1400, &
use_bnd1_cldtf_for_h2o_bands, &
calc_co2_tfs_monthly, &
calc_ch4_tfs_monthly, &
calc_n2o_tfs_monthly, &
calc_co2_tfs_on_first_step, &
use_current_co2_for_tf, &
co2_tf_calc_intrvl, &
co2_tf_time_displacement, &
calc_ch4_tfs_on_first_step, &
use_current_ch4_for_tf, &
ch4_tf_calc_intrvl, &
ch4_tf_time_displacement, &
calc_n2o_tfs_on_first_step, &
use_current_n2o_for_tf, &
n2o_tf_calc_intrvl, &
n2o_tf_time_displacement
!---------------------------------------------------------------------
!------- public data ------
!---------------------------------------------------------------------
!------- private data ------
!---------------------------------------------------------------------
! apcm, bpcm capphi coefficients for NBLY bands.
! atpcm, btpcm cappsi coefficients for NBLY bands.
! acomb random "a" parameter for NBLY bands.
! bcomb random "b" parameter for NBLY bands.
!---------------------------------------------------------------------
real, dimension (:), allocatable :: apcm, bpcm, atpcm, btpcm,&
acomb, bcomb
!-------------------------------------------------------------------
! the following longwave tables are retained for the life of the
! run.
!-------------------------------------------------------------------
type (longwave_tables3_type), save :: tabsr
type (longwave_tables1_type), save :: tab1, tab2, tab3, tab1w
type (longwave_tables2_type), save :: tab1a, tab2a, tab3a
!-------------------------------------------------------------------
!
!--------------------------------------------------------------------
integer, parameter :: no_combined_bands = 8
real, dimension(no_combined_bands) :: band_no_start, band_no_end
integer, dimension(NBLY_CKD-1) :: cld_indx_table
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
real, dimension (:), allocatable :: c1b7, c2b7
integer, dimension (:), allocatable :: cld_indx
!---------------------------------------------------------------------
! miscellaneous variables.
!---------------------------------------------------------------------
integer :: nbly ! number of frequency bands for exact
! cool-to-space computations.
integer :: nbtrge, nbtrg
integer :: ixprnlte
integer :: ks, ke
logical :: do_co2_tf_calc = .true.
logical :: do_ch4_tf_calc = .true.
logical :: do_n2o_tf_calc = .true.
logical :: do_co2_tf_calc_init = .true.
logical :: do_ch4_tf_calc_init = .true.
logical :: do_n2o_tf_calc_init = .true.
integer :: month_of_co2_tf_calc = 0
integer :: month_of_ch4_tf_calc = 0
integer :: month_of_n2o_tf_calc = 0
!----------------------------------------------------------------------
!----------------------------------------------------------------------
contains
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!
! PUBLIC SUBROUTINES
!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!#####################################################################
!
!
! Subroutine to initialize longwave radiation
!
!
! This subroutine initializes longwave radiation. It includes the
! prescribed gas band coefficients, initializes gas optical depth,
! longwave tables, and allocate cloud related variables in the
! longwave spectrum.
!
!
! call sealw99_init (latb, lonb, pref, Lw_tables)
!
!
! 2d array of model latitudes at cell corners [radians]
!
!
! 2d array of model longitudes at cell corners [radians]
!
!
! array containing two reference pressure profiles [pascals]
!
!
! lw_tables_type variable containing various longwave
! table specifiers needed by radiation_diag_mod.
!
!
!
subroutine sealw99_init (latb, lonb, pref, Lw_tables)
!---------------------------------------------------------------------
! sealw99_init is the constructor for sealw99_mod.
!---------------------------------------------------------------------
real, dimension(:,:), intent(in) :: latb, lonb
real, dimension(:,:), intent(in) :: pref
type(lw_table_type), intent(inout) :: Lw_tables
!---------------------------------------------------------------------
! intent(in) variables:
!
! lonb 2d array of model longitudes on cell corners
! [ radians ]
! latb 2d array of model latitudes at cell corners
! [ radians ]
! pref array containing two reference pressure profiles
! for use in defining transmission functions
! [ Pa ]
!
! intent(inout)
!
! Lw_tables lw_table_type variable which holds much of the
! relevant data used in the longwave radiation
! parameterization
!
!---------------------------------------------------------------------
!---------------------------------------------------------------------
! local variables:
real, dimension (no_combined_bands) :: band_no_start_rsb, &
band_no_start_ckd, &
band_no_end_rsb, &
band_no_end_ckd
data band_no_start_ckd / 1, 25, 41, 42, 43, 44, 46, 47 /
data band_no_end_ckd / 24,40, 41, 42, 43, 45, 46, 47 /
data band_no_start_rsb / 1, 5, 9, 10, 11, 12, 14, 15 /
data band_no_end_rsb / 4, 8, 9, 10, 11, 13, 14, 15 /
real, dimension (NBLY_CKD-1) :: cld_indx_table_lwclde, &
cld_indx_table_rsb
data cld_indx_table_lwclde /40*1, 2, 2, 2, 3, 4, 5, 6 /
data cld_indx_table_rsb / 47*1 /
real, dimension(NBLY_RSB) :: apcm_n, bpcm_n, &
atpcm_n, btpcm_n, &
acomb_n, bcomb_n
real, dimension(NBLY_CKD) :: apcm_c, bpcm_c, atpcm_c, &
btpcm_c, acomb_c, bcomb_c
real, dimension(size(pref,1) ) :: plm
real, dimension (NBCO215) :: cent, del
integer :: unit, ierr, io, k, n, nn, logunit
integer :: ioffset
real :: prnlte
integer :: kmax, kmin
integer :: inrad
real :: dum
character(len=4) :: gas_name
!---------------------------------------------------------------------
! local variables:
!
! band_no_start_rsb
! band_no_start_ckd
! band_no_end_rsb
! band_no_end_ckd
! cld_indx_table_lwclde
! cld_indx_table_rsb
! apcm_n
! bpcm_n
! atpcm_n
! btpcm_n
! acomb_n
! bcomb_n
! apcm_c
! bpcm_c
! atpcm_c
! btpcm_c
! acomb_c
! bcomb_c
! plm
! cent
! del
! unit
! ierr
! io
! k,n,m,nn
! ioffset
! prnlte
! kmax
! kmin
! inrad
! dum
!
!--------------------------------------------------------------------
!---------------------------------------------------------------------
! 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
!-----------------------------------------------------------------------
! read namelist.
#ifdef INTERNAL_FILE_NML
read (input_nml_file, nml=sealw99_nml, iostat=io)
ierr = check_nml_error(io,"sealw99_nml")
#else
!-----------------------------------------------------------------------
if ( file_exist('input.nml')) then
unit = open_namelist_file ( )
ierr=1; do while (ierr /= 0)
read (unit, nml=sealw99_nml, iostat=io, end=10)
ierr = check_nml_error(io,'sealw99_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=sealw99_nml)
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
kmax = size(pref,1) - 1 ! radiation grid size
ks = 1
ke = kmax
!---------------------------------------------------------------------
! be sure that the radiation time step has been defined before using
! it.
!---------------------------------------------------------------------
if (Rad_control%lw_rad_time_step_iz ) then
else
call error_mesg ('sealw99_mod', &
'must define lw_rad_time_step before using it', FATAL)
endif
!---------------------------------------------------------------------
! call check_tf_interval to verify that the namelist input variables
! related to co2 are consistent.
!---------------------------------------------------------------------
gas_name = 'co2 '
call check_tf_interval (gas_name, co2_tf_calc_intrvl, &
calc_co2_tfs_on_first_step, &
calc_co2_tfs_monthly, &
use_current_co2_for_tf)
!--------------------------------------------------------------------
! define the radiation_control_type variable components related to
! co2 tf calculation.
!--------------------------------------------------------------------
Rad_control%co2_tf_calc_intrvl = co2_tf_calc_intrvl
Rad_control%use_current_co2_for_tf = use_current_co2_for_tf
Rad_control%calc_co2_tfs_monthly = &
calc_co2_tfs_monthly
Rad_control%calc_co2_tfs_on_first_step = &
calc_co2_tfs_on_first_step
Rad_control%co2_tf_time_displacement = co2_tf_time_displacement
Rad_control%co2_tf_calc_intrvl_iz = .true.
Rad_control%use_current_co2_for_tf_iz = .true.
Rad_control%calc_co2_tfs_on_first_step_iz = .true.
Rad_control%calc_co2_tfs_monthly_iz = .true.
Rad_control%co2_tf_time_displacement_iz = .true.
!---------------------------------------------------------------------
! call check_tf_interval to verify that the namelist input variables
! related to ch4 are consistent.
!---------------------------------------------------------------------
gas_name = 'ch4 '
call check_tf_interval (gas_name, ch4_tf_calc_intrvl, &
calc_ch4_tfs_on_first_step, &
calc_ch4_tfs_monthly, &
use_current_ch4_for_tf)
!--------------------------------------------------------------------
! define the radiation_control_type variable components related to
! ch4 tf calculation.
!--------------------------------------------------------------------
Rad_control%ch4_tf_calc_intrvl = ch4_tf_calc_intrvl
Rad_control%use_current_ch4_for_tf = use_current_ch4_for_tf
Rad_control%calc_ch4_tfs_on_first_step = &
calc_ch4_tfs_on_first_step
Rad_control%calc_ch4_tfs_monthly = &
calc_ch4_tfs_monthly
Rad_control%ch4_tf_time_displacement = ch4_tf_time_displacement
Rad_control%ch4_tf_calc_intrvl_iz = .true.
Rad_control%use_current_ch4_for_tf_iz = .true.
Rad_control%calc_ch4_tfs_on_first_step_iz = .true.
Rad_control%calc_ch4_tfs_monthly_iz = .true.
Rad_control%ch4_tf_time_displacement_iz = .true.
!---------------------------------------------------------------------
! call check_tf_interval to verify that the namelist input variables
! related to n2o are consistent.
!---------------------------------------------------------------------
gas_name = 'n2o '
call check_tf_interval (gas_name, n2o_tf_calc_intrvl, &
calc_n2o_tfs_on_first_step, &
calc_n2o_tfs_monthly, &
use_current_n2o_for_tf)
!--------------------------------------------------------------------
! define the radiation_control_type variable components related to
! n2o tf calculation.
!--------------------------------------------------------------------
Rad_control%n2o_tf_calc_intrvl = n2o_tf_calc_intrvl
Rad_control%use_current_n2o_for_tf = use_current_n2o_for_tf
Rad_control%calc_n2o_tfs_on_first_step = &
calc_n2o_tfs_on_first_step
Rad_control%calc_n2o_tfs_monthly = &
calc_n2o_tfs_monthly
Rad_control%n2o_tf_time_displacement = n2o_tf_time_displacement
Rad_control%n2o_tf_calc_intrvl_iz = .true.
Rad_control%use_current_n2o_for_tf_iz = .true.
Rad_control%calc_n2o_tfs_on_first_step_iz = .true.
Rad_control%calc_n2o_tfs_monthly_iz = .true.
Rad_control%n2o_tf_time_displacement_iz = .true.
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
call longwave_params_init
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
if (trim(continuum_form) == 'ckd2.1' .or. &
trim(continuum_form) == 'ckd2.4' .or. &
trim(continuum_form) == 'mt_ckd1.0' .or. &
trim(continuum_form) == 'rsb' .or. &
trim(continuum_form) == 'none' ) then
else
call error_mesg ( 'sealw99_mod', &
'continuum_form is not specified correctly', FATAL)
endif
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
if (trim(linecatalog_form) == 'hitran_1992' .or. &
trim(linecatalog_form) == 'hitran_2000' ) then
else
call error_mesg ( 'sealw99_mod', &
'linecatalog_form is not specified correctly', FATAL)
endif
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
Lw_control%continuum_form = continuum_form
Lw_control%linecatalog_form = linecatalog_form
Lw_control%do_ch4lbltmpint = do_ch4lbltmpint
Lw_control%do_n2olbltmpint = do_n2olbltmpint
Lw_control%do_ch4lbltmpint_iz = .true.
Lw_control%do_n2olbltmpint_iz = .true.
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
if (trim(Lw_control%continuum_form) == 'ckd2.1' .or. &
trim(Lw_control%continuum_form) == 'ckd2.4' ) then
Lw_parameters%offset = 32
NBLY = NBLY_CKD
else if (trim(Lw_control%continuum_form) == 'rsb' ) then
Lw_parameters%offset = 0
NBLY = NBLY_RSB
endif
Lw_parameters%offset_iz = .true.
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
if (trim(Lw_control%linecatalog_form) == 'hitran_1992' ) then
if (trim(Lw_control%continuum_form) == 'ckd2.1' .or. &
trim(Lw_control%continuum_form) == 'ckd2.4' ) then
inrad = open_namelist_file ('INPUT/h2ocoeff_ckd_speccombwidebds_hi92')
read (inrad,9000) dum
read (inrad,9000) dum
read (inrad,9000) dum ! ckd capphi coeff for 560-800 band
read (inrad,9000) dum ! ckd cappsi coeff for 560-800 band
read (inrad,9000) dum ! ckd capphi coeff for 560-800 band
read (inrad,9000) dum ! ckd cappsi coeff for 560-800 band
read (inrad,9000) dum ! lo freq of 560-800 band
read (inrad,9000) dum ! hi freq of 560-800 band
! ckd rndm coeff for 40 bands (160-560) and 8 wide bands (560-1400)
read (inrad,9000) (acomb_c(k),k=1,NBLY_CKD)
read (inrad,9000) (bcomb_c(k),k=1,NBLY_CKD)
read (inrad,9000) (apcm_c(k),k=1,NBLY_CKD)
read (inrad,9000) (bpcm_c(k),k=1,NBLY_CKD)
read (inrad,9000) (atpcm_c(k),k=1,NBLY_CKD)
read (inrad,9000) (btpcm_c(k),k=1,NBLY_CKD)
else if (trim(Lw_control%continuum_form) == 'rsb' ) then
inrad = open_namelist_file ('INPUT/h2ocoeff_rsb_speccombwidebds_hi92')
read (inrad,9000) dum
read (inrad,9000) dum
read (inrad,9000) dum ! rsb capphi coeff for 560-800 band
read (inrad,9000) dum ! rsb cappsi coeff for 560-800 band
read (inrad,9000) dum ! rsb capphi coeff for 560-800 band
read (inrad,9000) dum ! rsb cappsi coeff for 560-800 band
read (inrad,9000) dum ! lo freq of 560-800 band
read (inrad,9000) dum ! hi freq of 560-800 band
read (inrad,9000) dum
! rsb rndm coeff for 8 comb bands (160-560) and 8 wide bands (560-1400)
read (inrad,9000) (acomb_n(k),k=1,NBLY_RSB)
read (inrad,9000) (bcomb_n(k),k=1,NBLY_RSB)
read (inrad,9000) (apcm_n(k),k=1,NBLY_RSB)
read (inrad,9000) (bpcm_n(k),k=1,NBLY_RSB)
read (inrad,9000) (atpcm_n(k),k=1,NBLY_RSB)
read (inrad,9000) (btpcm_n(k),k=1,NBLY_RSB)
endif
else if (trim(Lw_control%linecatalog_form) == 'hitran_2000' ) then
if (trim(Lw_control%continuum_form) == 'ckd2.1' .or. &
trim(Lw_control%continuum_form) == 'ckd2.4' ) then
inrad = open_namelist_file ('INPUT/h2ocoeff_ckd_speccombwidebds_hi00')
read (inrad,9000) dum
read (inrad,9000) dum
read (inrad,9000) dum ! ckd capphi coeff for 560-800 band
read (inrad,9000) dum ! ckd cappsi coeff for 560-800 band
read (inrad,9000) dum ! ckd capphi coeff for 560-800 band
read (inrad,9000) dum ! ckd cappsi coeff for 560-800 band
read (inrad,9000) dum ! lo freq of 560-800 band
read (inrad,9000) dum ! hi freq of 560-800 band
! ckd rndm coeff for 40 bands (160-560) and 8 wide bands (560-1400)
read (inrad,9000) (acomb_c(k),k=1,NBLY_CKD)
read (inrad,9000) (bcomb_c(k),k=1,NBLY_CKD)
read (inrad,9000) (apcm_c(k),k=1,NBLY_CKD)
read (inrad,9000) (bpcm_c(k),k=1,NBLY_CKD)
read (inrad,9000) (atpcm_c(k),k=1,NBLY_CKD)
read (inrad,9000) (btpcm_c(k),k=1,NBLY_CKD)
else if (trim(Lw_control%continuum_form) == 'rsb' ) then
inrad = open_namelist_file ('INPUT/h2ocoeff_rsb_speccombwidebds_hi00')
read (inrad,9000) dum
read (inrad,9000) dum
read (inrad,9000) dum ! rsb capphi coeff for 560-800 band
read (inrad,9000) dum ! rsb cappsi coeff for 560-800 band
read (inrad,9000) dum ! rsb capphi coeff for 560-800 band
read (inrad,9000) dum ! rsb cappsi coeff for 560-800 band
read (inrad,9000) dum ! lo freq of 560-800 band
read (inrad,9000) dum ! hi freq of 560-800 band
read (inrad,9000) dum
! rsb rndm coeff for 8 comb bands (160-560) and 8 wide bands (560-1400)
read (inrad,9000) (acomb_n(k),k=1,NBLY_RSB)
read (inrad,9000) (bcomb_n(k),k=1,NBLY_RSB)
read (inrad,9000) (apcm_n(k),k=1,NBLY_RSB)
read (inrad,9000) (bpcm_n(k),k=1,NBLY_RSB)
read (inrad,9000) (atpcm_n(k),k=1,NBLY_RSB)
read (inrad,9000) (btpcm_n(k),k=1,NBLY_RSB)
endif
endif
9000 format (5e14.6)
call close_file (inrad)
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
!---------------------------------------------------------------------
! define the number of separate bands in the 1200 - 1400 cm-1
! region. this region may be broken into 1, 2 , 4, 10 or 20 individ-
! ual bands, or it may remain part of the 0-160, 1200-2200 band.
!---------------------------------------------------------------------
nbtrge = no_h2o_bands_1200_1400
nbtrg = no_h2o_bands_1200_1400
Lw_parameters%NBTRG = nbtrg
Lw_parameters%NBTRGE = nbtrge
!---------------------------------------------------------------------
! set flag indicating these values have been initialized.
!---------------------------------------------------------------------
Lw_parameters%NBTRG_iz = .true.
Lw_parameters%NBTRGE_iz = .true.
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
call lw_gases_stdtf_init (pref)
call optical_path_init (pref)
call longwave_tables_init (Lw_tables, tabsr, &
tab1, tab2, tab3, tab1w, &
tab1a, tab2a, tab3a)
if (Lw_control%do_co2_iz) then
if (do_nlte .and. .not. Lw_control%do_co2) then
call error_mesg ('sealw99_mod', &
' cannot activate nlte when co2 not active as radiative gas',&
FATAL)
endif
else ! (do_co2_iz)
call error_mesg ('sealw99_mod', &
'do_co2 not yet defined', FATAL)
endif ! (do_co2_iz)
!---------------------------------------------------------------------
! define pressure-dependent index values used in the infrared
! radiation code. by this manner, the coefficients are defined
! at execution time (not dependent on the choice of vertical
! layers)
! prnlte : pressure (mb) below which non-LTE code (Nlte.F) affects
! CO2 transmissivities
!---------------------------------------------------------------------
if (do_nlte) then
prnlte = 0.1
!--------------------------------------------------------------------
! abort execution if trying to run with modified radiation grid
! code must be added to properly map plm (on the model grid)
! to the radiation grid. if simply dropping upper levels, then is
! fairly easy. abort here so that problem may be addressed and to
! prevent uncertified results.
! solution is likely to be to pass in plm on radiation grid, whatever
! that is.
!--------------------------------------------------------------------
kmin = 1
plm (kmin) = 0.
do k=kmin+1,kmax
plm (k) = 0.5*(pref (k-1,1) + pref (k,1))
end do
plm (kmax+1) = pref (kmax+1,1)
!--------------------------------------------------------------------
! convert pressure specification for bottom (flux) pressure level
! for nlte calculation into an index (ixprnlte)
!!! CAN THE MODEL TOP BE AT A PRESSURE THAT IS NOT ZERO ??
!!! if not, then must define plm(ks) always to be 0.0
!!! implications for lw_gas tf calcs ??
! kmax = size(plm,1) - 1
!-------------------------------------------------------------------
ixprnlte = 1
do k=ks+1, kmax
if (plm(k)*1.0E-02 .LT. prnlte) then
ixprnlte = k-ks+1
else
exit
endif
end do
!---------------------------------------------------------------------
! allocate and obtain elements of the source function for bands in
! the 15 um range (used in nlte)
!---------------------------------------------------------------------
ioffset = Lw_parameters%offset
allocate ( c1b7 (NBCO215) )
allocate ( c2b7 (NBCO215) )
do n=1,NBCO215
cent(n) = 0.5E+00*(Lw_tables%bdlocm(n+8+ioffset) + &
Lw_tables%bdhicm(n+8+ioffset))
del (n) = Lw_tables%bdhicm(n+8+ioffset) - &
Lw_tables%bdlocm(n+8+ioffset)
c1b7(n) = (3.7412E-05)*cent(n)*cent(n)*cent(n)*del(n)
c2b7(n) = (1.4387E+00)*cent(n)
end do
endif
!--------------------------------------------------------------------
!
!--------------------------------------------------------------------
call gas_tf_init (pref)
!--------------------------------------------------------------------
!
!--------------------------------------------------------------------
call longwave_clouds_init
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
if (trim(Lw_control%continuum_form) == 'ckd2.1' .or. &
trim(Lw_control%continuum_form) == 'ckd2.4' ) then
band_no_start = band_no_start_ckd
band_no_end = band_no_end_ckd
NBLY = NBLY_CKD
else if (trim(Lw_control%continuum_form) == 'rsb' ) then
band_no_start = band_no_start_rsb
band_no_end = band_no_end_rsb
NBLY = NBLY_RSB
endif
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
Lw_parameters%NBLY = NBLY
Lw_parameters%NBLY_iz = .true.
Lw_control%do_lwcldemiss = do_lwcldemiss
Lw_control%do_lwcldemiss_iz = .true.
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
if (do_lwcldemiss) then
cld_indx_table = cld_indx_table_lwclde
Cldrad_control%nlwcldb = 7
else
cld_indx_table = cld_indx_table_rsb
Cldrad_control%nlwcldb = 1
endif
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
allocate (apcm(NBLY))
allocate (bpcm(NBLY))
allocate (atpcm(NBLY))
allocate (btpcm(NBLY))
allocate (acomb(NBLY))
allocate (bcomb(NBLY))
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
if (trim(Lw_control%continuum_form) == 'ckd2.1' .or. &
trim(Lw_control%continuum_form) == 'ckd2.4' ) then
apcm = apcm_c
atpcm = atpcm_c
bpcm = bpcm_c
btpcm = btpcm_c
acomb = acomb_c
bcomb = bcomb_c
else if (trim(Lw_control%continuum_form) == 'rsb' ) then
apcm = apcm_n
atpcm = atpcm_n
bpcm = bpcm_n
btpcm = btpcm_n
acomb = acomb_n
bcomb = bcomb_n
endif
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
allocate (cld_indx(6+NBTRGE) )
do nn=1,6+NBTRGE
if (nn > Cldrad_control%nlwcldb) then
cld_indx(nn) = Cldrad_control%nlwcldb
else
cld_indx(nn) = nn
endif
end do
if (NBTRGE == 0 .and. .not. use_bnd1_cldtf_for_h2o_bands) then
call error_mesg ('sealw99_mod', &
'must use band1 cld tfs for the 1200-1400 cm(-1) bands when &
& they are included in the 1200-2200 cm(-1) band', FATAL)
endif
if (NBTRGE > 0 .and. use_bnd1_cldtf_for_h2o_bands) then
cld_indx(7:) = 1
endif
!--------------------------------------------------------------------
!
!---------------------------------------------------------------------
call longwave_fluxes_init
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
module_is_initialized = .true.
!---------------------------------------------------------------------
end subroutine sealw99_init
!####################################################################
subroutine sealw99_time_vary (Rad_time, Rad_gases)
type(time_type), intent(in) :: Rad_time
type(radiative_gases_type), intent(inout) :: Rad_gases
logical :: calc_co2, calc_n2o, calc_ch4
character(len=4) :: gas_name
integer :: year, month, day, hour, minute, &
second
!---------------------------------------------------------------------------
call get_control_gas_tf (calc_co2, calc_ch4, calc_n2o)
!----------------------------------------------------------------------
!
!--------------------------------------------------------------------
call lw_gases_stdtf_time_vary
!----------------------------------------------------------------------
!
!--------------------------------------------------------------------
if (Rad_gases%time_varying_ch4 .or. &
Rad_gases%time_varying_n2o) then
if (Rad_gases%time_varying_ch4 .and. .not. calc_ch4) then
call error_mesg ('sealw99_mod', &
' if ch4 amount is to vary in time, ch4 tfs must be '//&
'recalculated', FATAL)
endif
if (Rad_gases%time_varying_n2o .and. .not. calc_n2o) then
call error_mesg ('sealw99_mod', &
' if n2o amount is to vary in time, n2o tfs must be '//&
'recalculated', FATAL)
endif
endif
!----------------------------------------------------------------------
!
!--------------------------------------------------------------------
if (Rad_gases%time_varying_co2) then
if (.not. calc_co2) then
call error_mesg ('sealw99_mod', &
' if co2 amount is to vary in time, co2 tfs must be '//&
'recalculated', FATAL)
endif
endif
!----------------------------------------------------------------------
! if ch4 is activated in this job, varying in time, and
! calculation of ch4 tfs are requested, call obtain_gas_tfs to
! define the tfs.
!--------------------------------------------------------------------
if (Rad_gases%time_varying_ch4) then
if (Lw_control%do_ch4) then
if (do_ch4_tf_calc) then
gas_name = 'ch4 '
call obtain_gas_tfs (gas_name, Rad_time, &
Rad_gases%Ch4_time, &
ch4_tf_calc_intrvl,&
Rad_gases%ch4_tf_offset, &
calc_ch4_tfs_on_first_step, &
calc_ch4_tfs_monthly, &
month_of_ch4_tf_calc, &
Rad_gases%ch4_for_next_tf_calc, &
Rad_gases%ch4_for_last_tf_calc, &
do_ch4_tf_calc, do_ch4_tf_calc_init)
endif ! (do_ch4_tf_calc)
endif ! (do_ch4)
!---------------------------------------------------------------------
! if ch4 is not time-varying and it is the initial call to sealw99,
! call ch4_time_vary to calculate the tfs. set flags to indicate
! the calculation has been done.
!---------------------------------------------------------------------
else ! (time_varying_ch4)
if (Lw_control%do_ch4 .and. do_ch4_tf_calc ) then
call ch4_time_vary (Rad_gases%rrvch4)
do_ch4_tf_calc = .false.
do_ch4_tf_calc_init = .false.
else if (.not. Lw_control%do_ch4) then
do_ch4_tf_calc = .false.
do_ch4_tf_calc_init = .false.
endif
endif ! (time_varying_ch4)
!----------------------------------------------------------------------
! if n2o is activated in this job, varying in time, and
! calculation of n2o tfs are requested, call obtain_gas_tfs to
! define the tfs.
!--------------------------------------------------------------------
if (Rad_gases%time_varying_n2o) then
if (Lw_control%do_n2o) then
if (do_n2o_tf_calc) then
gas_name = 'n2o '
call obtain_gas_tfs (gas_name, Rad_time, &
Rad_gases%N2o_time, &
n2o_tf_calc_intrvl,&
Rad_gases%n2o_tf_offset, &
calc_n2o_tfs_on_first_step, &
calc_n2o_tfs_monthly, &
month_of_n2o_tf_calc, &
Rad_gases%n2o_for_next_tf_calc, &
Rad_gases%n2o_for_last_tf_calc, &
do_n2o_tf_calc, do_n2o_tf_calc_init)
endif ! (do_n2o_tf_calc)
endif ! (do_n2o)
!---------------------------------------------------------------------
! if n2o is not time-varying and it is the initial call to sealw99,
! call n2o_time_vary to calculate the tfs. set flags to indicate
! the calculation has been done.
!---------------------------------------------------------------------
else
if (Lw_control%do_n2o .and. do_n2o_tf_calc) then
call n2o_time_vary (Rad_gases%rrvn2o)
do_n2o_tf_calc = .false.
do_n2o_tf_calc_init = .false.
else if (.not. Lw_control%do_n2o) then
do_n2o_tf_calc = .false.
do_n2o_tf_calc_init = .false.
endif
endif ! (time_varying_n2o)
!----------------------------------------------------------------------
! if co2 is activated in this job, varying in time, and
! calculation of co2 tfs are requested, call obtain_gas_tfs to
! define the tfs.
!--------------------------------------------------------------------
if (Rad_gases%time_varying_co2) then
if (Lw_control%do_co2) then
if (do_co2_tf_calc) then
gas_name = 'co2 '
call obtain_gas_tfs (gas_name, Rad_time, &
Rad_gases%Co2_time, &
co2_tf_calc_intrvl,&
Rad_gases%co2_tf_offset, &
calc_co2_tfs_on_first_step, &
calc_co2_tfs_monthly, &
month_of_co2_tf_calc, &
Rad_gases%co2_for_next_tf_calc, &
Rad_gases%co2_for_last_tf_calc, &
do_co2_tf_calc, do_co2_tf_calc_init)
endif ! (do_co2_tf_calc)
endif ! (do_co2)
!---------------------------------------------------------------------
! if co2 is not time-varying and it is the initial call to sealw99,
! call co2_time_vary to calculate the tfs. set flags to indicate
! the calculation has been done.
!---------------------------------------------------------------------
else
! interactive co2 mod for radiation calculation
! here it's hardcoded to recompute co2 TF on the 1st of each month
if (Rad_gases%use_model_supplied_co2) then
call get_date (Rad_time, year, month, day, hour, minute,&
second)
if (day == 1 .and. hour == 0 .and. minute == 0 .and. &
second == 0) then
call co2_time_vary (Rad_gases%rrvco2)
Rad_gases%co2_for_last_tf_calc = Rad_gases%rrvco2
do_co2_tf_calc_init = .false.
else
if (do_co2_tf_calc_init) then
call co2_time_vary (Rad_gases%co2_for_last_tf_calc)
do_co2_tf_calc_init = .false.
endif
endif
else !(Rad_gases%use_model_supplied_co2)
if (Lw_control%do_co2 .and. do_co2_tf_calc) then
call co2_time_vary (Rad_gases%rrvco2)
do_co2_tf_calc = .false.
do_co2_tf_calc_init = .false.
else if (.not. Lw_control%do_co2) then
do_co2_tf_calc = .false.
do_co2_tf_calc_init = .false.
endif
endif !(Rad_gases%use_model_supplied_co2)
endif ! (time_varying_co2)
!----------------------------------------------------------------------
!
!--------------------------------------------------------------------
if ((Lw_control%do_co2 .and. calc_co2) .or. &
(Lw_control%do_ch4 .and. calc_ch4) .or. &
(Lw_control%do_n2o .and. calc_n2o)) then
call lw_gases_stdtf_dealloc
endif
!------------------------------------------------------------------------
end subroutine sealw99_time_vary
!#####################################################################
subroutine sealw99_endts (Rad_gases_tv)
type(radiative_gases_type), intent(in) :: Rad_gases_tv
if (Rad_gases_tv%time_varying_ch4) then
if (Lw_control%do_ch4) then
if (.not. calc_ch4_tfs_on_first_step) then
do_ch4_tf_calc = .true.
endif
endif
endif
if (Rad_gases_tv%time_varying_n2o) then
if (Lw_control%do_n2o) then
if (.not. calc_n2o_tfs_on_first_step) then
do_n2o_tf_calc = .true.
endif
endif
endif
if (Rad_gases_tv%time_varying_co2) then
if (Lw_control%do_co2) then
if (.not. calc_co2_tfs_on_first_step) then
do_co2_tf_calc = .true.
endif
endif
endif
end subroutine sealw99_endts
!#####################################################################
!#####################################################################
!
!
! Subroutine to calculate longwave radiation flux and heating rate.
!
!
! This subroutine calculates longwave radiation flux and heating rate
! based on the simplified exchange method. It also provides diagnostics
! for the longwave radiation and cloud.
!
!
! call sealw99 (is, ie, js, je, Atmos_input, Rad_gases, &
! Aerosol, Aerosol_props, Cldrad_props, Cld_spec, &
! Lw_output, Lw_diagnostics)
!
!
! starting subdomain i indice of data in the physics_window being
! integrated
!
!
! ending subdomain i indice of data in the physics_window being
! integrated
!
!
! starting subdomain j indice of data in the physics_window being
! integrated
!
!
! ending subdomain j indice of data in the physics_window being
! integrated
!
!
! atmos_input_type variable containing the atmospheric
! input fields needed by the radiation package
!
!
! radiative_gases_type variable containing the radi-
! ative gas input fields needed by the radiation
! package
!
!
! Aerosol climatological input data to longwave radiation
!
!
! Aerosol radiative properties
!
!
! cldrad_properties_type variable containing the
! cloud radiative property input fields needed by the
! radiation package
!
!
! Cloud specification type contains cloud microphysical, geometrical,
! and distribution properties in a model column.
!
!
! lw_output_type variable containing longwave
! radiation output data
!
!
! lw_diagnostics_type variable containing diagnostic
! longwave output used by the radiation diagnostics
! module
!
!
!
subroutine sealw99 (is, ie, js, je, Rad_time, Atmos_input, Rad_gases, &
Aerosol, Aerosol_props, Cldrad_props, Cld_spec, &
Aerosol_diags, Lw_output, Lw_diagnostics, &
including_aerosols)
!---------------------------------------------------------------------
! sealw99 is the longwave driver subroutine.
!
! references:
!
! (1) schwarzkopf, m. d., and s. b. fels, "the simplified
! exchange method revisited: an accurate, rapid method for
! computation of infrared cooling rates and fluxes," journal
! of geophysical research, 96 (1991), 9075-9096.
!
! (2) schwarzkopf, m. d., and s. b. fels, "improvements to the
! algorithm for computing co2 transmissivities and cooling
! rates," journal geophysical research, 90 (1985) 10541-10550.
!
! (3) fels, s.b., "simple strategies for inclusion of voigt
! effects in infrared cooling calculations," application
! optics, 18 (1979), 2634-2637.
!
! (4) fels, s. b., and m. d. schwarzkopf, "the simplified exchange
! approximation: a new method for radiative transfer
! calculations," journal atmospheric science, 32 (1975),
! 1475-1488.
!
! author: m. d. schwarzkopf
!
! revised: 10/7/93
!
! certified: radiation version 1.0
!---------------------------------------------------------------------
integer, intent(in) :: is, ie, js, je
type(time_type), intent(in) :: Rad_time
type(atmos_input_type), intent(in) :: Atmos_input
type(radiative_gases_type), intent(inout) :: Rad_gases
type(aerosol_type), intent(in) :: Aerosol
type(aerosol_properties_type), intent(inout) :: Aerosol_props
type(aerosol_diagnostics_type), intent(inout) :: Aerosol_diags
type(cldrad_properties_type), intent(in) :: Cldrad_props
type(cld_specification_type), intent(in) :: Cld_spec
type(lw_output_type), intent(inout) :: Lw_output
type(lw_diagnostics_type), intent(inout) :: Lw_diagnostics
logical, intent(in) :: including_aerosols
!-----------------------------------------------------------------------
! intent(in) variables:
!
! is,js,ie,je
! Atmos_input
! Rad_gases
! Aerosol
! Cldrad_props
! Cld_spec
!
! intent(inout) variables:
!
! Aerosol_props
! Lw_output
! Lw_diagnostics
!
!---------------------------------------------------------------------
!---------------------------------------------------------------------
! local variables:
real, dimension (size(Atmos_input%press,1), &
size(Atmos_input%press,2), &
size(Atmos_input%press,3), &
Cldrad_control%nlwcldb) :: &
cldtf
real, dimension (size(Atmos_input%press,1), &
size(Atmos_input%press,2), &
size(Atmos_input%press,3)) :: &
cnttaub1, cnttaub2, cnttaub3, co21c, &
co21r, heatem, overod, tmp1, to3cnt
real, dimension (size(Atmos_input%press,1), &
size(Atmos_input%press,2), 2) :: &
emspec
real, dimension (size(Atmos_input%press,1), &
size(Atmos_input%press,2)) :: &
co21c_KEp1, co21r_KEp1
real, dimension (size(Atmos_input%press,1), &
size(Atmos_input%press,2), &
size(Atmos_input%press,3), 3) :: &
contdg
real, dimension (size(Atmos_input%press,1), &
size(Atmos_input%press,2), &
size(Atmos_input%press,3) -1) :: &
pdflux
real, dimension (size(Atmos_input%press,1), &
size(Atmos_input%press,2)) :: &
flx1e1cf, gxctscf
real, dimension (size(Atmos_input%press,1), &
size(Atmos_input%press,2), &
size(Atmos_input%press,3)) :: &
e1ctw1, e1ctw2, &
emisdg, flxcf, &
heatemcf, flx, to3dg, tcfc8
real, dimension (size(Atmos_input%press,1), &
size(Atmos_input%press,2), &
size(Atmos_input%press,3) - 1) :: &
cts_sum, cts_sumcf
real, dimension (size(Atmos_input%press,1), &
size(Atmos_input%press,2), &
size(Atmos_input%press,3), NBTRGE) :: &
emisdgf, tch4n2oe
real, dimension (size(Atmos_input%press,1), &
size(Atmos_input%press,2), &
size(Atmos_input%press,3), &
8):: &
sorc
real, dimension (size(Atmos_input%press,1), &
size(Atmos_input%press,2), &
size(Atmos_input%press,3), &
6+NBTRGE ) :: &
source_band, dsrcdp_band
real, dimension (size(Atmos_input%press,1), &
size(Atmos_input%press,2), &
size(Atmos_input%press,3) , &
6+NBTRGE ) :: &
trans_band1, trans_band2
real, dimension (size(Atmos_input%press,1), &
size(Atmos_input%press,2), &
6+NBTRGE ) :: &
trans_b2d1, trans_b2d2
real, dimension (size(Atmos_input%press,1), &
size(Atmos_input%press,2), 2, NBTRGE) :: &
emspecf
real, dimension (size(Atmos_input%press,1), &
size(Atmos_input%press,2), &
size(Atmos_input%press,3)-1 ) :: &
pdfinv, heatra_save
real, dimension (size(Atmos_input%press,1), &
size(Atmos_input%press,2), &
size(Atmos_input%press,3) ) :: flxnet_save
type(lw_clouds_type) :: Lw_clouds
type(optical_path_type) :: Optical
type(gas_tf_type) :: Gas_tf
integer :: ix, jx, kx
integer :: k, kp, m, j
integer :: kk, i, l
integer :: nprofiles, nnn
!---------------------------------------------------------------------
! local variables:
!
! cldtf
! cnttaub1
! cnttaub2
! cnttaub3
! co21c transmission function for the 560-800 cm-1 band
! from levels k through KE+1 to level k. used for flux
! at level k arising from exchange with levels k
! through KE+1. includes co2 (from tables), h2o (after
! multiplication with over).
! co21diag transmission function for co2 only in the 560-800
! cm-1 band, from levels k to k, where k ranges from
! KS to KE+1.
! co21r transmission function for the 560-800 cm-1 band
! from level k to levels k+1 through KE+1. used for
! flux at levels k+1 through KE+1 arising from
! exchange with level k. includes co2 (from tables),
! h2o (after multiplication with over).
! dsorc15
! dsorc93
! dsorcb1
! dsorcb2
! dsorcb3
! dt4
! emiss
! heatem
! overod
! sorc15 planck function for 560-800 cm-1 bands (sum over
! bands 9 and 10).
! t4
! tmp1 temporary array, used for computational purposes
! in various places. should have no consequences
! beyond the immediate module wherein it is defined.
! tmp2 temporary arrays, similar to tmp1
! to3cnt transmission function for the 990-1070 cm-1 band
! including o3(9.6 um) + h2o continuum (no lines)
! and possibly cfcs.
! ch41c
! n2o1c
! n2o17c
! emspec
! s1a
! flxge1
! co21c_KEp1
! co21r_KEp1
! contdg
! cfc_tf
! pdflux
! flx1e1cf
! flxge1cf
! gxctscf
! emissb
! e1cts1
! e1cts2
! e1ctw1
! e1ctw2
! soe2
! soe3
! soe4
! soe5
! emisdg
! flxcf
! heatemcf
! flx
! to3dg
! taero8
! taero8kp
! totaer_tmp
! tcfc8
! cts_sum
! cts_sumcf
! emissbf
! e1cts1f
! e1cts2f
! emisdgf
! emissf
! tch4n2oe
! flx1e1fcf
! flxge1f
! flxge1fcf
! sorc planck function, at model temperatures, for all
! bands; used in cool-to-space calculations.
! source_band
! dsrcdp_band
! trans_band1
! trans_band2
! trans_b2d1
! trans_b2d2
! emspecf
! pdfinv
! Lw_clouds
! Optical
! Gas_tf
! calc_co2
! calc_n2o
! calc_ch4
! ch4_vmr
! n2o_vmr
! co2_vmr
! ix,jx,kx
! n,k,kp,m,j
!
!---------------------------------------------------------------------
!---------------------------------------------------------------------
! be sure module is initialized.
!---------------------------------------------------------------------
if (.not. module_is_initialized ) then
call error_mesg( 'sealw99_mod', &
'module has not been initialized', FATAL )
endif
!----------------------------------------------------------------------
!
!--------------------------------------------------------------------
kx = size(Atmos_input%press,3) - 1
ix = ie-is+1
jx = je-js+1
!--------------------------------------------------------------------
! call sealw99_alloc to allocate component arrays of
! lw_diagnostics_type variables.
!----------------------------------------------------------------------
if ( .not. associated (Lw_diagnostics%flx1e1)) then
call sealw99_alloc (ix, jx, kx, Lw_diagnostics)
else
Lw_diagnostics%flx1e1 = 0.
Lw_diagnostics%cts_out = 0.
Lw_diagnostics%cts_outcf = 0.
Lw_diagnostics%gxcts = 0.
Lw_diagnostics%excts = 0.
Lw_diagnostics%exctsn = 0.
Lw_diagnostics%fctsg = 0.
Lw_diagnostics%fluxn = 0.
if (Rad_control%do_totcld_forcing) then
Lw_diagnostics%fluxncf = 0.
endif
Lw_diagnostics%flx1e1f = 0.
endif
!----------------------------------------------------------------------
!
!--------------------------------------------------------------------
call optical_path_setup (is, ie, js, je, Atmos_input, Rad_gases, &
Aerosol, Aerosol_props, Aerosol_diags, &
Optical, including_aerosols)
!--------------------------------------------------------------------
! call co2coef to compute some co2 temperature and pressure
! interpolation quantities and to compute temperature-corrected co2
! transmission functions (co2spnb and co2nbl).
!-------------------------------------------------------------------
call co2coef (Atmos_input, Gas_tf)
!----------------------------------------------------------------------
! call co2_source_calc to calculate the source function.
!-----------------------------------------------------------------------
call co2_source_calc (Atmos_input, Rad_gases, sorc, Gas_tf, &
source_band, dsrcdp_band)
if (Cldrad_control%do_ica_calcs) then
nprofiles = Cldrad_control%nlwcldb
heatra_save = 0.
flxnet_save = 0.0
else
nprofiles = 1
endif
do nnn = 1, nprofiles
! reinitialize these values (done in sealw99_alloc)
Lw_diagnostics%flx1e1 = 0.
Lw_diagnostics%cts_out = 0.
Lw_diagnostics%cts_outcf = 0.
Lw_diagnostics%gxcts = 0.
Lw_diagnostics%excts = 0.
Lw_diagnostics%exctsn = 0.
Lw_diagnostics%fctsg = 0.
Lw_diagnostics%fluxn (:,:,:,:) = 0.0
if (Rad_control%do_totcld_forcing) then
Lw_diagnostics%fluxncf(:,:,:,:) = 0.0
endif
if (NBTRGE > 0) then
Lw_diagnostics%flx1e1f = 0.
end if
!----------------------------------------------------------------------
! call cldtau to compute cloud layer transmission functions for
! all layers.
!-------------------------------------------------------------------
call cldtau (nnn, Cldrad_props, Cld_spec, Lw_clouds)
!---------------------------------------------------------------------
! BEGIN LEVEL KS CALCULATIONS
!---------------------------------------------------------------------
!-----------------------------------------------------------------
! compute co2 560-800 cm-1, ch4 and n2o 1200-1400 cm-1 trans-
! mission functions and n2o 560-670 cm-1 transmission functions
! appropriate for level KS.
!---------------------------------------------------------------------
call transcolrow (Gas_tf, KS, KS, KS, KE+1, KS+1, KE+1, &
co21c, co21r, tch4n2oe)
!---------------------------------------------------------------------
! go into optical_path_mod to obtain the optical path functions
! needed for use from level KS.
! to3cnt contains values in the 990 - 1070 cm-1 range (ozone, water
! vapor continuum, aerosol, cfc).
! overod contains values in the 560 - 800 cm-1 range (water vapor
! lines, continuum, aerosol, 17 um n2o band, cfc).
! cnttaub1 is continuum band 4 (water vapor continuum, aerosol, cfc)
! cnttaub2 is continuum band 5 (water vapor continuum, aerosol, cfc)
! cnttaub3 is continuum band 7 (water vapor continuum, aerosol, cfc)
! the 15um band transmission functions between levels KS and KS+1
! are stored in overod and co2nbl; they will not be overwritten,
! as long as calculations are made for pressure levels increasing
! from KS.
!---------------------------------------------------------------------
call optical_trans_funct_from_KS (Gas_tf, to3cnt, overod, &
Optical, cnttaub1, cnttaub2, &
cnttaub3, including_aerosols)
!-----------------------------------------------------------------------
! compute cloud transmission functions between level KS and all
! other levels.
!-----------------------------------------------------------------------
call cloud (KS, Cldrad_props, Cld_spec, Lw_clouds, cldtf)
!-----------------------------------------------------------------------
! obtain exact cool-to-space for water and co2, and approximate
! cool-to-space for co2 and o3.
!----------------------------------------------------------------------
call cool_to_space_exact (cldtf, Atmos_input, Optical, Gas_tf, &
sorc, to3cnt, Lw_diagnostics, cts_sum, &
cts_sumcf, gxctscf, including_aerosols)
!----------------------------------------------------------------------
! compute the emissivity fluxes for k=KS.
!! trans_band1:
! index 1 = e1flx
! index 2 = co21r*overod
! index 3 = cnttaub1
! index 4 = cnttaub2
! index 5 = to3cnt
! index 6 = cnttaub3
! index 7 = e1flxf
!! trans_band2:
! index 1 = emiss
! index 2 = co21c*overod
! index 3 = cnttaub1
! index 4 = cnttaub2
! index 5 = to3cnt
! index 6 = cnttaub3
! index 7 = emissf
!----------------------------------------------------------------------
call e1e290 (Atmos_input, e1ctw1, e1ctw2, trans_band1, &
trans_band2, Optical, tch4n2oe, &
source_band(:,:,:,1), Lw_diagnostics, cldtf, &
cld_indx, flx1e1cf, tcfc8, including_aerosols)
do kk = KS,KE+1
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kk,2) = co21r (i,j,kk)* &
overod(i,j,kk)
end do
end do
end do
do kk = KS,KE+1
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kk,3) = cnttaub1(i,j,kk)
end do
end do
end do
do kk = KS,KE+1
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kk,4) = cnttaub2(i,j,kk)
end do
end do
end do
do kk = KS,KE+1
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kk,5) = to3cnt (i,j,kk)
end do
end do
end do
do kk = KS,KE+1
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kk,6) = cnttaub3(i,j,kk)
end do
end do
end do
!! trans_band2:
! index 1 = emiss
! index 2 = co21c*overod
! index 3 = cnttaub1
! index 4 = cnttaub2
! index 5 = to3cnt
! index 6 = cnttaub3
! index 7 = emissf
do kk = KS+1,KE+1
do j = 1,size(trans_band2(:,:,:,:),2)
do i = 1,size(trans_band2(:,:,:,:),1)
trans_band2(i,j,kk,2) = co21c(i,j,kk)* &
overod(i,j,kk)
end do
end do
end do
do l = 3,6
do kk = KS+1,KE+1
do j = 1,size(trans_band2(:,:,:,:),2)
do i = 1,size(trans_band2(:,:,:,:),1)
trans_band2(i,j,kk,l) = trans_band1(i,j,kk,l)
end do
end do
end do
end do
!----------------------------------------------------------------------
! the following is a rewrite of the original code largely to
! eliminate three-dimensional arrays. the code works on the
! following principles. let k be a fixed flux level and kp be
! a varying flux level, then
!
! flux(k) = sum(deltab(kp)*tau(kp,k)) for kp=KS,KE+1.
!
! if we assume a symmetrical array tau(k,kp)=tau(kp,k), we can
! break down the calculations for k=KS,KE+1 as follows:
!
! flux(k) = sum(deltab(kp)*tau(kp,k)) for kp=k+1,KE+1 (1)
!
! flux(kp) = (deltab(k )*tau(kp,k)) for kp=k+1,KE+1. (2)
!
! plus deltab(k)*tau(k,k) for all k.
!
! if we compute a one-dimensional array tauod(kp) for
! kp=k+1,KE+1, equations (1) and (2) become:
!
! tauod(kp) = tau(kp,k) (3)
!
! flux (k ) = sum(deltab(kp)*tauod(kp)) for kp=k+1,KE+1 (4)
!
! flux (kp) = (deltab(k )*tauod(kp)) for kp=k+1,KE+1 (5)
!
! where tau(k,k) and nearby layer terms are handled separately.
!
! compute fluxes at level k = KS
! compute the terms for flux at levels KS+1 to KE+1 from level KS.
! compute terms for flux at level KS from level KS.
! compute the terms for flux at level KS due to levels KP from KS+1
! to KE+1.
!
!-----------------------------------------------------------------------
call longwave_fluxes_ks (source_band, trans_band1, dsrcdp_band, &
trans_band2, cldtf, cld_indx , &
Lw_diagnostics)
do kk = KS,KE+1
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kk,1) = e1ctw2(i,j,kk)
end do
end do
end do
!-----------------------------------------------------------------------
! compute approximate cool-to-space heating rates for 1 wide band
! in the 15um range (560-800 cm-1) (ctsco2) and for 1 band in
! the 9.6 um band (ctso3).
!----------------------------------------------------------------------
call cool_to_space_approx ( Atmos_input%pflux, source_band, &
trans_band1, cldtf, cld_indx , &
Lw_diagnostics, e1ctw1)
!-----------------------------------------------------------------------
! perform flux calculations for the flux levels KS+1 to KE-1. calcu-
! lations for flux levels KE and KE+1 are done separately, as all
! calculations are special cases or nearby layers.
!----------------------------------------------------------------------
do k=KS+1,KE-1
!--------------------------------------------------------------------
! compute co2 560-800 cm-1, ch4 and n2o 1200-1400 cm-1 trans-
! mission functions and n2o 560-670 cm-1 transmission functions
! appropriate for level k.
!--------------------------------------------------------------------
call transcolrow (Gas_tf, k, k, k, KE+1, k+1, KE+1, &
co21c, co21r, tch4n2oe)
!-------------------------------------------------------------------
! the 15 um band transmission functions between levels k and k+1
! are stored in overod and co2nbl; they will not be overwritten,
! as long as calculations are made for pressure levels increasing
! from k.
!---------------------------------------------------------------------
call optical_trans_funct_k_down (Gas_tf, k, to3cnt, overod, &
Optical, including_aerosols)
!-----------------------------------------------------------------------
! compute cloud transmission functions between level k and all
! other levels greater or equal to k.
!---------------------------------------------------------------------
call cloud (k,Cldrad_props,Cld_spec, Lw_clouds, cldtf)
!-----------------------------------------------------------------------
! compute the exchange terms in the flux equation (except the
! nearby layer (k,k) terms, done later).
!! trans_band1:
! index 1 = emissb
! index 2 = co21r*overod
! index 3 = contodb1
! index 4 = contodb2
! index 5 = to3cnt
! index 6 = contodb3
! index 7 = emissbf
!! trans_band2:
! index 1 = emiss
! index 2 = co21c*overod
! index 3 = contodb1
! index 4 = contodb2
! index 5 = to3cnt
! index 6 = contodb3
! index 7 = emissf
!---------------------------------------------------------------------
call e290 (Atmos_input, k, trans_band2, trans_band1, &
Optical, tch4n2oe, tcfc8, including_aerosols)
do kp=k,KE
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kp+1,3) = cnttaub1(i,j,kp+1 )/cnttaub1(i,j,k )
end do
end do
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kp+1,4) = cnttaub2(i,j,kp+1 )/cnttaub2(i,j,k )
end do
end do
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kp+1,6) = cnttaub3(i,j,kp+1 )/cnttaub3(i,j,k )
end do
end do
end do
do kk = k+1,KE+1
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kk,2) = co21r(i,j,kk)* &
overod(i,j,kk)
end do
end do
end do
do kk = k+1,KE+1
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kk,5) = to3cnt(i,j,kk)
end do
end do
end do
do kk = k+1,KE+1
do j = 1,size(trans_band2(:,:,:,:),2)
do i = 1,size(trans_band2(:,:,:,:),1)
trans_band2(i,j,kk,2) = co21c(i,j,kk)* &
overod(i,j,kk)
end do
end do
end do
do l = 3,6
do kk = k+1,KE+1
do j = 1,size(trans_band2(:,:,:,:),2)
do i = 1,size(trans_band2(:,:,:,:),1)
trans_band2(i,j,kk,l) = trans_band1(i,j,kk,l)
end do
end do
end do
end do
!-----------------------------------------------------------------------
! compute the terms for flux at levels k+1 to KE+1 from level k.
! compute the terms for flux at level k due to levels
! kp from k+1 to KE+1.
!----------------------------------------------------------------------
call longwave_fluxes_k_down (k, dsrcdp_band, trans_band1, &
trans_band2, cldtf, cld_indx, &
Lw_diagnostics )
end do ! (end of k=KS+1,KE-1 loop)
!-----------------------------------------------------------------------
! compute remaining flux terms. these include:
! 1) the (k,k) terms, for pressure levels k from KS+1 to KE-1
! (the KS,KS term was handled earlier);
! 2) terms for pressure level KE. these include the (KE,KE) term,
! computed as in (1), and the (KE,KE+1) and (KE+1,KE) terms,
! computed somewhat differently from the similar terms at
! higher levels, owing to the proximity to the surface layer
! KE+1;
! 3) the term for pressure level KE+1 (the (KE+1,KE+1 term).
!
! compute k=KE case. since the kp loop is length one, many
! simplifications occur. the co2 quantities and the emissivity
! quantities) are computed in the nbl section. therefore, we want
! to compute over, to3cnt, and contod; according to our notation
! over(:,:,KE), to3cnt(:,:,KE), and contod(:,:,KE). the boundary
! layer and nearby layer corrections to the transmission functions
! are obtained above. the following ratios are used in various nbl
! nbl calculations. the remaining calculations are for:
!
! 1) the (k,k) terms, k=KS+1,KE-1;
! 2) the (KE,KE ) term;
! 3) the (KE,KE+1 ) term;
! 4) the (KE+1,KE ) term;
! 5) the (KE+1,KE+1) term.
!
! each is uniquely handled. different flux terms are computed
! differently the fourth section obtains water transmission
! functions used in q(approximate) calculations and also makes nbl
! corrections:
!
! 1) emiss (:,:) is the transmission function matrix obtained
! using E2spec;
!
! 2) "nearby layer" corrections (emiss(i,i)) are obtained
! using E3v88;
!
! 3) special values at the surface (emiss(KE,KE+1),
! emiss(KE+1,KE), emiss(KE+1,KE+1)) are calculated.
!
!
! compute temperature and/or scaled amount indices and residuals
! for nearby layer and special layer lookup tables.
!
! calculation for special cases (KE,KE+1) and (KE+1,KE)
!
! compute co2 560-800 cm-1, ch4 and n2o 1200-1400 cm-1 trans-
! mission functions and n2o 560-670 cm-1 transmission functions
! appropriate for level KE. if activated, save ch4n2o tf term.
!-------------------------------------------------------------------
call transcolrow (Gas_tf, KE, KE, KE, KE+1, KE+1, KE+1, &
co21c, co21r, tch4n2oe)
!----------------------------------------------------------------------
! get optical path terms for KE
!----------------------------------------------------------------------
call optical_trans_funct_KE (Gas_tf, to3cnt, Optical, overod, &
including_aerosols)
do j = 1,size(trans_b2d1(:,:,:),2)
do i = 1,size(trans_b2d1(:,:,:),1)
trans_b2d1(i,j,3 ) = cnttaub1(i,j,KE+1)/cnttaub1(i,j,KE )
end do
end do
do j = 1,size(trans_b2d1(:,:,:),2)
do i = 1,size(trans_b2d1(:,:,:),1)
trans_b2d1(i,j,4 ) = cnttaub2(i,j,KE+1)/cnttaub2(i,j,KE )
end do
end do
do j = 1,size(trans_b2d1(:,:,:),2)
do i = 1,size(trans_b2d1(:,:,:),1)
trans_b2d1(i,j,6 ) = cnttaub3(i,j,KE+1)/cnttaub3(i,j,KE )
end do
end do
!-----------------------------------------------------------------------
! compute cloud transmission functions between level KE and KE and
! KE+1
!----------------------------------------------------------------------
call cloud (KE, Cldrad_props, Cld_spec, Lw_clouds, cldtf)
!--------------------------------------------------------------------
! compute mean temperature in the "nearby layer" between a flux
! level and the first data level below the flux level (tpl1) or the
! first data level above the flux level (tpl2)
!---------------------------------------------------------------------
call esfc (Atmos_input, emspec, Optical, emspecf, &
tch4n2oe, tcfc8)
!----------------------------------------------------------------------
! compute nearby layer transmission functions for 15 um band, cont-
! inuum bands, and 9.3 um band in subroutine Nearbylyrtf. trans-
! mission functions for the special cases (KE,KE+1) and (KE+1,KE)
! are also computed for the 15 um band.
!! trans_band1:
! index 1 = emspec(KS+1)
! index 2 = co21c_KEp1
! index 3 = contodb1
! index 4 = contodb2
! index 5 = to3cnt
! index 6 = contodb3
! index 7 = emspecf(KS+1)
!! trans_band2:
! index 1 = emspec(KS)
! index 2 = co21r_KEp1
! index 3 = contodb1
! index 4 = contodb2
! index 5 = to3cnt
! index 6 = contodb3
! index 7 = emspecf(KS)
!----------------------------------------------------------------------
call trans_sfc (Gas_tf, Atmos_input, overod, co21c_KEp1, &
co21r_KEp1)
do j = 1,size(trans_b2d1(:,:,:),2)
do i = 1,size(trans_b2d1(:,:,:),1)
trans_b2d1(i,j,1) = emspec(i,j,KS+1)
end do
end do
do j = 1,size(trans_b2d1(:,:,:),2)
do i = 1,size(trans_b2d1(:,:,:),1)
trans_b2d1(i,j,2) = co21c_KEp1(i,j)
end do
end do
do j = 1,size(trans_b2d1(:,:,:),2)
do i = 1,size(trans_b2d1(:,:,:),1)
trans_b2d1(i,j,5) = to3cnt(i,j,KE+1)
end do
end do
do m=1,NBTRGE
do j = 1,size(trans_b2d1(:,:,:),2)
do i = 1,size(trans_b2d1(:,:,:),1)
trans_b2d1(i,j,6+m) = emspecf(i,j,KS+1,m)
end do
end do
end do
do j = 1,size(trans_b2d2(:,:,:),2)
do i = 1,size(trans_b2d2(:,:,:),1)
trans_b2d2(i,j,1) = emspec(i,j,KS)
end do
end do
do j = 1,size(trans_b2d2(:,:,:),2)
do i = 1,size(trans_b2d2(:,:,:),1)
trans_b2d2(i,j,2) = co21r_KEP1(i,j)
end do
end do
do kk = 3,6
do j = 1,size(trans_b2d2(:,:,:),2)
do i = 1,size(trans_b2d2(:,:,:),1)
trans_b2d2(i,j,kk) = trans_b2d1(i,j,kk)
end do
end do
end do
do m=1,NBTRGE
do j = 1,size(trans_b2d2(:,:,:),2)
do i = 1,size(trans_b2d2(:,:,:),1)
trans_b2d2(i,j,6+m) = emspecf(i,j,KS,m)
end do
end do
end do
!-----------------------------------------------------------------------
! obtain fluxes for the two terms (KE,KE+1) and (KE+1,KE), both
! using the same cloud transmission functions (from layer KE)
!----------------------------------------------------------------------
call longwave_fluxes_KE_KEp1 (dsrcdp_band, trans_b2d1, &
trans_b2d2, cldtf, cld_indx, &
Lw_diagnostics )
!---------------------------------------------------------------------
! call enear to calculate emissivity arrays
!----------------------------------------------------------------------
call enear (Atmos_input, emisdg, Optical, emisdgf , tch4n2oe, &
tcfc8)
!-------------------------------------------------------------------
! obtain optical path transmission functions for diagonal terms
!------------------------------------------------------------------
call optical_trans_funct_diag (Atmos_input, contdg, to3dg, &
Optical)
!-----------------------------------------------------------------------
! compute cloud transmission functions between level KE+1 and KE+1
!------------------------------------------------------------------
call cloud (KE+1, Cldrad_props, Cld_spec, Lw_clouds, cldtf)
!----------------------------------------------------------------------
! compute nearby layer transmission functions for 15 um band, cont-
! inuum bands, and 9.3 um band in subroutine Nearbylyrtf. trans-
! mission functions for the special cases (KE,KE+1) and (KE+1,KE)
! are also computed for the 15 um band.
!----------------------------------------------------------------------
call trans_nearby (Gas_tf, Atmos_input, overod, co21c)
do kk = ks+1,KE+1
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kk,1) = emisdg(i,j,kk)
end do
end do
end do
do kk = ks+1,KE+1
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kk,2) = co21c(i,j,kk)
end do
end do
end do
do kk = ks+1,KE+1
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kk,3) = contdg(i,j,kk,1)
end do
end do
end do
do kk = ks+1,KE+1
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kk,4) = contdg(i,j,kk,2)
end do
end do
end do
do kk = ks+1,KE+1
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kk,5) = to3dg(i,j,kk)
end do
end do
end do
do kk = ks+1,KE+1
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kk,6) = contdg(i,j,kk,3)
end do
end do
end do
do m=1,NBTRGE
do kk = ks+1,KE+1
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kk,6+m) = emisdgf(i,j,kk,m)
end do
end do
end do
end do
!-----------------------------------------------------------------------
! obtain fluxes for the diagonal terms at all levels.
!-----------------------------------------------------------------------
call longwave_fluxes_diag (dsrcdp_band, trans_band1, cldtf , &
cld_indx, Lw_diagnostics )
!--------------------------------------------------------------------
! sum up fluxes over bands
!-----------------------------------------------------------------------
if (Rad_control%do_totcld_forcing) then
call longwave_fluxes_sum (is, ie, js, je, flx, NBTRGE, &
Lw_diagnostics, flxcf)
else
call longwave_fluxes_sum (is, ie, js, je, flx, NBTRGE, &
Lw_diagnostics)
endif
Lw_output%bdy_flx(:,:,1) = Lw_output%bdy_flx(:,:,1) + &
Lw_diagnostics%fluxn(:,:,1,3) + &
Lw_diagnostics%fluxn(:,:,1,4) + &
Lw_diagnostics%fluxn(:,:,1,5) + &
Lw_diagnostics%fluxn(:,:,1,6)
Lw_output%bdy_flx(:,:,2) = Lw_output%bdy_flx(:,:,2) + &
Lw_diagnostics%fluxn(:,:,1,4)
Lw_output%bdy_flx(:,:,3) = Lw_output%bdy_flx(:,:,3) + &
Lw_diagnostics%fluxn(:,:,ke+1,3) + &
Lw_diagnostics%fluxn(:,:,ke+1,4) + &
Lw_diagnostics%fluxn(:,:,ke+1,5) + &
Lw_diagnostics%fluxn(:,:,ke+1,6)
Lw_output%bdy_flx(:,:,4) = Lw_output%bdy_flx(:,:,4) + &
Lw_diagnostics%fluxn(:,:,ke+1,4)
! Lw_output%bdy_flx = 1.0E-03*Lw_output%bdy_flx
if (nnn == 1) then ! need do only once
if (Rad_control%do_totcld_forcing) then
Lw_output%bdy_flx_clr(:,:,1) = Lw_diagnostics%fluxncf(:,:,1,3) + &
Lw_diagnostics%fluxncf(:,:,1,4) + &
Lw_diagnostics%fluxncf(:,:,1,5) + &
Lw_diagnostics%fluxncf(:,:,1,6)
Lw_output%bdy_flx_clr(:,:,2) = Lw_diagnostics%fluxncf(:,:,1,4)
Lw_output%bdy_flx_clr(:,:,3) = Lw_diagnostics%fluxncf(:,:,ke+1,3) + &
Lw_diagnostics%fluxncf(:,:,ke+1,4) + &
Lw_diagnostics%fluxncf(:,:,ke+1,5) + &
Lw_diagnostics%fluxncf(:,:,ke+1,6)
Lw_output%bdy_flx_clr(:,:,4) = Lw_diagnostics%fluxncf(:,:,ke+1,4)
Lw_output%bdy_flx_clr = 1.0E-03*Lw_output%bdy_flx_clr
endif
endif
!-----------------------------------------------------------------------
! compute emissivity heating rates.
!-----------------------------------------------------------------------
do kk = ks,ke
do j = 1,size(pdfinv(:,:,:),2)
do i = 1,size(pdfinv(:,:,:),1)
pdfinv(i,j,kk) = 1.0/(Atmos_input%pflux(i,j,kk+ks+1-(ks)) - &
Atmos_input%pflux(i,j,kk))
end do
end do
end do
do kk = KS,KE
do j = 1,size(heatem(:,:,:),2)
do i = 1,size(heatem(:,:,:),1)
heatem(i,j,kk) = (radcon_mks*(flx(i,j,kk+KS+1-(KS)) - &
flx(i,j,kk))*pdfinv(i,j,kk))*1.0e-03
end do
end do
end do
if (Rad_control%do_totcld_forcing) then
do kk = KS,KE
do j = 1,size(heatemcf(:,:,:),2)
do i = 1,size(heatemcf(:,:,:),1)
heatemcf(i,j,kk) = (radcon_mks*(flxcf(i,j,kk+KS+1-(KS)) - &
flxcf(i,j,kk))*pdfinv(i,j,kk)*1.0e-03)
end do
end do
end do
endif
!-----------------------------------------------------------------------
! compute total heating rates.
!-----------------------------------------------------------------------
!--------------------------------------------------------------------
! cts_sum is the sum of the values from cool_to_space_exact and
! the values defined here in cool_to_space_approx. it will be used
! by longwave_driver_mod.
!--------------------------------------------------------------------
do kk = 1,size(cts_sum(:,:,:),3)
do j = 1,size(cts_sum(:,:,:),2)
do i = 1,size(cts_sum(:,:,:),1)
cts_sum(i,j,kk) = ((((((cts_sum(i,j,kk) - &
Lw_diagnostics%cts_out(i,j,kk,2)) - &
Lw_diagnostics%cts_out(i,j,kk,5) )- &
Lw_diagnostics%cts_out(i,j,kk,1)) - &
Lw_diagnostics%cts_out(i,j,kk,3)) - &
Lw_diagnostics%cts_out(i,j,kk,4))- &
Lw_diagnostics%cts_out(i,j,kk,6))
end do
end do
end do
if (Rad_control%do_totcld_forcing) then
do kk = 1,size(cts_sumcf(:,:,:),3)
do j = 1,size(cts_sumcf(:,:,:),2)
do i = 1,size(cts_sumcf(:,:,:),1)
cts_sumcf(i,j,kk) = ((((((cts_sumcf(i,j,kk) - &
Lw_diagnostics%cts_outcf(i,j,kk,2)) - &
Lw_diagnostics%cts_outcf(i,j,kk,5)) - &
Lw_diagnostics%cts_outcf(i,j,kk,1)) - &
Lw_diagnostics%cts_outcf(i,j,kk,3)) - &
Lw_diagnostics%cts_outcf(i,j,kk,4) ) - &
Lw_diagnostics%cts_outcf(i,j,kk,6))
end do
end do
end do
endif
do kk = KS,KE
do j = 1,size(Lw_output%heatra(:,:,:),2)
do i = 1,size(Lw_output%heatra(:,:,:),1)
Lw_output%heatra(i,j,kk) = heatem(i,j,kk) + &
cts_sum (i,j,kk)
end do
end do
end do
if (nnn == 1) then ! only need to do once
if (Rad_control%do_totcld_forcing) then
do kk = KS,KE
do j = 1,size(Lw_output%heatracf(:,:,:),2)
do i = 1,size(Lw_output%heatracf(:,:,:),1)
Lw_output%heatracf(i,j,kk) = heatemcf(i,j,kk) + &
cts_sumcf(i,j,kk)
end do
end do
end do
endif
endif ! (nnn == 1)
!-----------------------------------------------------------------------
! compute the flux at each flux level using the flux at the
! top (flx1e1 + gxcts) and the integral of the heating rates.
!---------------------------------------------------------------------
do kk = KS,KE
do j = 1,size(pdflux(:,:,:),2)
do i = 1,size(pdflux(:,:,:),1)
pdflux(i,j,kk) = Atmos_input%pflux(i,j,kk+KS+1-(KS)) - &
Atmos_input%pflux(i,j,kk)
end do
end do
end do
do j = 1,size(Lw_output%flxnet(:,:,:),2)
do i = 1,size(Lw_output%flxnet(:,:,:),1)
Lw_output%flxnet(i,j,KS ) = Lw_diagnostics%flx1e1(i,j) + Lw_diagnostics%gxcts(i,j)
end do
end do
!---------------------------------------------------------------------
! convert values to mks (1.0e-03 factor)
!---------------------------------------------------------------------
do j = 1,size(Lw_diagnostics%gxcts(:,:),2)
do i = 1,size(Lw_diagnostics%gxcts(:,:),1)
Lw_diagnostics%gxcts(i,j) = 1.0e-03*Lw_diagnostics%gxcts(i,j)
end do
end do
do j = 1,size(Lw_diagnostics%flx1e1(:,:),2)
do i = 1,size(Lw_diagnostics%flx1e1(:,:),1)
Lw_diagnostics%flx1e1(i,j) = 1.0e-03*Lw_diagnostics%flx1e1(i,j)
end do
end do
if (NBTRGE> 0) then
do kk = 1,size(Lw_diagnostics%flx1e1f(:,:,:),3)
do j = 1,size(Lw_diagnostics%flx1e1f(:,:,:),2)
do i = 1,size(Lw_diagnostics%flx1e1f(:,:,:),1)
Lw_diagnostics%flx1e1f(i,j,kk) = &
1.0e-03*Lw_diagnostics%flx1e1f(i,j,kk)
end do
end do
end do
endif
!---------------------------------------------------------------------
! convert mks values to cgs (1.0e03 factor) so can be summed with
! cgs value.
!---------------------------------------------------------------------
do kk = KS,KE
do j = 1,size(tmp1(:,:,:),2)
do i = 1,size(tmp1(:,:,:),1)
tmp1(i,j,kk) = 1.0e03*Lw_output%heatra(i,j,kk)*pdflux(i,j,kk)/radcon_mks
end do
end do
end do
do k=KS+1,KE+1
do j = 1,size(Lw_output%flxnet(:,:,:),2)
do i = 1,size(Lw_output%flxnet(:,:,:),1)
Lw_output%flxnet(i,j,k) = Lw_output%flxnet(i,j,k-1) + tmp1(i,j,k-1)
end do
end do
enddo
if (nnn == 1) then ! only need to do once
if (Rad_control%do_totcld_forcing) then
do j = 1,size(Lw_output%flxnetcf(:,:,:),2)
do i = 1,size(Lw_output%flxnetcf(:,:,:),1)
Lw_output%flxnetcf(i,j,KS ) = flx1e1cf(i,j) + gxctscf(i,j)
end do
end do
!---------------------------------------------------------------------
! convert mks values to cgs (1.0e03 factor) so can be summed
! with cgs value.
!---------------------------------------------------------------------
do kk = KS,KE
do j = 1,size(tmp1(:,:,:),2)
do i = 1,size(tmp1(:,:,:),1)
tmp1(i,j,kk) = 1.0e03*Lw_output%heatracf(i,j,kk)*pdflux(i,j,kk)/radcon_mks
end do
end do
end do
do k=KS+1,KE+1
do j = 1,size(Lw_output%flxnetcf(:,:,:),2)
do i = 1,size(Lw_output%flxnetcf(:,:,:),1)
Lw_output%flxnetcf(i,j,k) = Lw_output%flxnetcf(i,j,k-1) + tmp1(i,j,k-1)
end do
end do
enddo
endif
endif ! (nnn == 1)
if (Cldrad_control%do_ica_calcs) then
heatra_save = heatra_save + Lw_output%heatra
flxnet_save = flxnet_save + Lw_output%flxnet
endif
end do ! (profiles loop)
if (Cldrad_control%do_ica_calcs) then
Lw_output%heatra = heatra_save / Float(nprofiles)
Lw_output%flxnet = flxnet_save / Float(nprofiles)
Lw_output%bdy_flx = 1.0e-03*Lw_output%bdy_flx/Float(nprofiles)
endif
!-----------------------------------------------------------------------
! call thickcld to perform "pseudo-convective adjustment" for
! maximally overlapped clouds, if desired.
!-----------------------------------------------------------------------
if (do_thick) then
call thickcld (Atmos_input%pflux, Cldrad_props, Cld_spec, &
Lw_output)
endif ! (do_thick)
!--------------------------------------------------------------------
! convert lw fluxes to mks units.
!---------------------------------------------------------------------
do kk = 1,size(Lw_output%flxnet(:,:,:),3)
do j = 1,size(Lw_output%flxnet(:,:,:),2)
do i = 1,size(Lw_output%flxnet(:,:,:),1)
Lw_output%flxnet(i,j,kk) = 1.0E-03*Lw_output%flxnet(i,j,kk)
end do
end do
end do
if (Rad_control%do_totcld_forcing) then
do kk = 1,size(Lw_output%flxnetcf(:,:,:),3)
do j = 1,size(Lw_output%flxnetcf(:,:,:),2)
do i = 1,size(Lw_output%flxnetcf(:,:,:),1)
Lw_output%flxnetcf(i,j,kk) = 1.0E-03*Lw_output%flxnetcf(i,j,kk)
end do
end do
end do
endif
!--------------------------------------------------------------------
! call lw_clouds_dealloc to deallocate component arrays of Lw_clouds.
!--------------------------------------------------------------------
call lw_clouds_dealloc (Lw_clouds)
!--------------------------------------------------------------------
! call gas_tf_dealloc to deallocate component arrays of Gas_tf.
!--------------------------------------------------------------------
call gas_tf_dealloc (Gas_tf)
!--------------------------------------------------------------------
! call optical_dealloc to deallocate component arrays of Optical.
!--------------------------------------------------------------------
call optical_dealloc (Optical, including_aerosols)
!--------------------------------------------------------------------
end subroutine sealw99
!#####################################################################
!
!
! sealw99_end is the destructor for sealw99_mod.
!
!
! sealw99_end is the destructor for sealw99_mod.
!
!
! call sealw99_end
!
!
!
subroutine sealw99_end
!---------------------------------------------------------------------
! sealw99_end is the destructor for sealw99_mod.
!---------------------------------------------------------------------
!---------------------------------------------------------------------
! be sure module is initialized.
!---------------------------------------------------------------------
if (.not. module_is_initialized ) then
call error_mesg( 'sealw99_mod', &
'module has not been initialized', FATAL )
endif
!---------------------------------------------------------------------
! call the destructor routines for the modules initialized by
! sealw99.
!---------------------------------------------------------------------
call gas_tf_end
call optical_path_end
call lw_gases_stdtf_end
call longwave_clouds_end
call longwave_fluxes_end
call longwave_tables_end
call longwave_params_end
!---------------------------------------------------------------------
! mark the module as uninitialized.
!---------------------------------------------------------------------
module_is_initialized = .false.
!---------------------------------------------------------------------
end subroutine sealw99_end
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!
! PRIVATE SUBROUTINES
!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!#####################################################################
!
!
! check_tf_interval verifies that requested tf calculation intervals
! are compatible with radiation time step
!
!
!
! check_tf_interval verifies the following relationships:
! 1) that the tf calculation interval is no smaller than the
! radiation time step;
! 2) that the tf calculation interval is an integral multiple of
! the radiation time step;
! 3) that the specification for calculating tfs on the first step
! of the job is done properly.
!
!
! call check_tf_interval (gas, gas_tf_calc_intrvl, &
! calc_gas_tfs_on_first_step, &
! use_current_gas_for_tf)
! call sealw99_alloc (ix, jx, kx, Lw_diagnostics)
!
!
! name associated with the gas
!
!
! time interval between recalculating transmission fumctions [ hours ]
!
!
! flag indicating if tfs are to be calculated only on first step
! of job
!
!
! flag indicating if gas mixing ratio at current time is to be used
! for calculation of gas tfs
!
!
!
!
!NOTE: THIS IS A PRIVATE SUBROUTINE.
!
subroutine check_tf_interval (gas, gas_tf_calc_intrvl, &
calc_gas_tfs_on_first_step, &
calc_gas_tfs_monthly, &
use_current_gas_for_tf)
!--------------------------------------------------------------------
character(len=4), intent(in) :: gas
real, intent(in) :: gas_tf_calc_intrvl
logical, intent(in) :: calc_gas_tfs_on_first_step, &
calc_gas_tfs_monthly, &
use_current_gas_for_tf
!
!---------------------------------------------------------------------
! if tfs are not being calculated on the first step, the requested
! gas transmission function recalculation interval must be greater
! than the radiation time step.
!---------------------------------------------------------------------
if (.not. calc_gas_tfs_on_first_step .and. &
.not. calc_gas_tfs_monthly) then
if (INT(3600.0*gas_tf_calc_intrvl) < &
Rad_control%lw_rad_time_step) then
call error_mesg ('sealw99_mod', &
trim(gas)// ' tf calculation interval must be greater&
& than or equal to the radiation time step', FATAL)
endif
!---------------------------------------------------------------------
! be sure that the tf calculation interval is an integral multiple
! of the radiation timestep.
!---------------------------------------------------------------------
if (mod(INT(3600.0*gas_tf_calc_intrvl), &
Rad_control%lw_rad_time_step) /= 0) then
call error_mesg ('sealw99_mod', &
trim(gas)//' transmission function calculation interval &
&must be integral multiple of radiation time step', FATAL)
endif
endif ! (.not. calc_gas_tfs_on_first_step)
!---------------------------------------------------------------------
! to calculate the tfs using the gas value at the start of the run,
! one must set use_current_gas_for_tf to .false, and set the
! gas_tf_time_displacement to 0.0, rather than setting
! use_current_gas_for_tf to .true.
!---------------------------------------------------------------------
if (calc_gas_tfs_on_first_step) then
if (use_current_gas_for_tf) then
call error_mesg ('sealw99_mod', &
'cannot specify use of current '//trim(gas)//' value&
& for tfs when calculating tfs on first step; instead &
&set use_current_'//trim(gas)//'_for_tf to false and set &
& '//trim(gas)//'_tf_time_displacement = 0.0', FATAL)
endif
endif
if (calc_gas_tfs_on_first_step .and. &
calc_gas_tfs_monthly) then
call error_mesg ( 'sealw99_mod', &
'cannot request calc of tfs both on first step and monthly',&
FATAL)
endif
!---------------------------------------------------------------------
end subroutine check_tf_interval
!####################################################################
!
!
! obtain_gas_tfs obtains the transmission functions for the requested
! gas
!
!
!
! obtain_gas_tfs performs the following functions:
! a) if time variation of the gas has begun at the current time:
! 1) defines how long the gas has been varying and whether
! the tfs are due to be recalculated at the current time;
! 2) if the tfs are not to be always recalculated on the first
! step:
! a) if this is a recalculation step;
! 1) call the routine to calculate the tfs for the input
! gas;
! 2) redefine the value of the gas mixing ratio used fro the
! last tf calculation to be the one just used;
! 3) set the flag indicating the need to initially calculate
! the tfs to .false.
! b) if this is not a recalculation step:
! 1) if this is the initial step, call the routine to calc-
! ulate the tfs for the input gas;
! 2) set the flag indicating the need to initially calculate
! the tfs to .false.
! 3) if the tfs are to be always calculated on the first step:
! a) call the routine to calculate the tfs for the input gas;
! b) redefine the value of the gas mixing ratio used from the
! last tf calculation to be the one just used;
! c) set the flag indicating the need to initially calculate
! the tfs to .false.
! b) if time variation of the gas has not begun at the current time:
! 1) if this is the initial call of the job, call the routine
! to calculate the tfs for the input gas;
! 2) set a flag to indicate that the initial call has been made.
!
!
! call obtain_gas_tfs (gas, Rad_time, Gas_time, gas_tf_calc_intrvl,&
! gas_tf_offset, calc_gas_tfs_on_first_step, &
! gas_for_next_tf_calc, gas_for_last_tf_calc, &
! do_gas_tf_calc, do_gas_tf_calc_init)
!
!
! name associated with the gas
!
!
! current model time [ time_type ]
!
!
! time since time variation of gas began [ time_type ]
!
!
! time interval between recalculating transmission fumctions [ hours ]
!
!
! time difference between current time and the time for which the
! tfs are calculated
!
!
! flag indicating if tfs are to be calculated only on first step
! of job
!
!
! value of gas mixing ratio to be used when tfs are next calculated
! [ no. / no. ]
!
!
! value of gas mixing ratio to be used when tfs were last calculated
! [ no. / no. ]
!
!
! if true, calculate gas tfs when alarm again goes off
!
!
! this variable is true initially to force calculation of the tfs on
! the first call of the job; it is then set to false
!
!
!
!NOTE: THIS IS A PRIVATE SUBROUTINE.
!
subroutine obtain_gas_tfs (gas, Rad_time, Gas_time, gas_tf_calc_intrvl,&
gas_tf_offset, calc_gas_tfs_on_first_step, &
calc_gas_tfs_monthly,month_of_gas_tf_calc, &
gas_for_next_tf_calc, gas_for_last_tf_calc, &
do_gas_tf_calc, do_gas_tf_calc_init)
!----------------------------------------------------------------------
character(len=4), intent(in) :: gas
type(time_type), intent(in) :: Rad_time, Gas_time
real, intent(in) :: gas_tf_calc_intrvl, &
gas_tf_offset
logical, intent(in) :: calc_gas_tfs_on_first_step
integer, intent(inout) :: month_of_gas_tf_calc
logical, intent(in) :: calc_gas_tfs_monthly
real, intent(inout) :: gas_for_next_tf_calc, &
gas_for_last_tf_calc
logical, intent(inout) :: do_gas_tf_calc,&
do_gas_tf_calc_init
!---------------------------------------------------------------------
! local variables:
type(time_type) :: Time_since_gas_start
integer :: seconds, days, minutes_from_start, alarm
integer :: year, month, day, hour, minute, second
character(len=4) :: chvers, chvers2, chvers3, chvers4, &
chvers5, chvers6
character(len=12) :: chvers10
character(len=16) :: chvers7
!---------------------------------------------------------------------
! local variables:
!
! Time_since_gas_start length of time that gas has been varying
! [ time_type ]
! seconds seconds component of Time_since_gas_start
! [ seconds ]
! days days component of Time_since_gas_start
! [ days ]
! minutes_from_start minutes since time variation started
! [ minutes ]
! alarm if alarm = 0, it is time to calculate the
! tfs [ minutes ]
! year year component of current time
! month month component of current time
! day day component of current time
! hour hour component of current time
! minute minute component of current time
! second second component of current time
! chvers, chversx characters used to output model variables
! through the error_mesg interface
!
!---------------------------------------------------------------------
!--------------------------------------------------------------------
! if gas variation is underway, define how long it has been since
! that started, and whether the current time is an integral multiple
! of the calculation interval from that gas starting time. put
! the current time into character variables for output via the
! error_mesg interface.
!--------------------------------------------------------------------
if (Rad_time >= Gas_time) then
Time_since_gas_start = Rad_time - Gas_time
call get_time (Time_since_gas_start, seconds, days)
call get_date (Rad_time, year, month, day, hour, minute, second)
write (chvers, '(i4)') year
write (chvers2, '(i4)') month
write (chvers3, '(i4)') day
write (chvers4, '(i4)') hour
write (chvers5, '(i4)') minute
write (chvers6, '(i4)') second
write (chvers10, '( f9.3)') gas_tf_offset
!---------------------------------------------------------------------
! if tfs are not automatically calculated on the first step of the
! job and if the current time is a desired recalculation time, call
! xxx_time_vary to do the calculation.
!---------------------------------------------------------------------
! if (.not. calc_gas_tfs_on_first_step) then
if (.not. calc_gas_tfs_on_first_step .and. &
.not. calc_gas_tfs_monthly) then
minutes_from_start = INT(days*1440.0 + real(seconds)/60.)
if (gas_tf_calc_intrvl /= 0.0) then
alarm = MOD (minutes_from_start, &
INT(gas_tf_calc_intrvl*60.0))
endif
if (alarm == 0) then
if (trim(gas) == 'ch4') then
call ch4_time_vary (gas_for_next_tf_calc)
write (chvers7, '(4pe15.7)') gas_for_next_tf_calc
else if (trim(gas) == 'n2o') then
call n2o_time_vary (gas_for_next_tf_calc)
write (chvers7, '(3pe15.7)') gas_for_next_tf_calc
else if (trim(gas) == 'co2') then
call co2_time_vary (gas_for_next_tf_calc)
write (chvers7, '(3pe15.7)') gas_for_next_tf_calc
endif
!--------------------------------------------------------------------
! redefine the value for the gas mixing ratio used fro the last tf
! calculation.
!--------------------------------------------------------------------
gas_for_last_tf_calc = gas_for_next_tf_calc
!---------------------------------------------------------------------
! if a record of the tf calculation path is desired, print out the
! relevant data.
!---------------------------------------------------------------------
if (verbose >= 1) then
if (gas_tf_offset /= 0.0) then
call error_mesg ('sealw99_mod', &
'calculating '//trim(gas)//' transmission functions&
& at time '//chvers//chvers2//chvers3//chvers4// &
chvers5// chvers6// ', using '//trim(gas)//' &
&mixing ratio of:' // chvers7 //', which&
& is the value '//chvers10// 'hours from current &
& time.', NOTE)
else
call error_mesg ('sealw99_mod', &
'calculating '//trim(gas)//' transmission functions&
& at time ' //chvers//chvers2//chvers3//chvers4//&
chvers5// chvers6// ', using '//trim(gas)//' &
&mixing ratio of:' // chvers7 //', which&
& is the value at the current time', NOTE)
endif
endif ! (verbose)
!---------------------------------------------------------------------
! set the flag to indicate that the initial tf calculation has been
! completed.
!---------------------------------------------------------------------
do_gas_tf_calc_init = .false.
!----------------------------------------------------------------------
! if alarm is not 0 and the tfs have not yet been calculated, call
! xxx_time_vary to do the calculation. set the flag appropriately.
!---------------------------------------------------------------------
else
if (do_gas_tf_calc_init) then
if (trim(gas) == 'ch4') then
call ch4_time_vary (gas_for_last_tf_calc)
write (chvers7, '(4pe15.7)') gas_for_last_tf_calc
else if (trim(gas) == 'n2o') then
call n2o_time_vary (gas_for_last_tf_calc)
write (chvers7, '(3pe15.7)') gas_for_last_tf_calc
else if (trim(gas) == 'co2') then
call co2_time_vary (gas_for_last_tf_calc)
write (chvers7, '(3pe15.7)') gas_for_last_tf_calc
endif
do_gas_tf_calc_init = .false.
!---------------------------------------------------------------------
! if a record of the tf calculation path is desired, print out the
! relevant data.
!---------------------------------------------------------------------
if (verbose >= 1) then
if (gas_tf_offset /= 0.0) then
call error_mesg ('sealw99_mod', &
'initial '//gas//' transmission function &
&calculation uses '//trim(gas)//' mixing ratio &
&of:'//chvers7//'.', NOTE)
else
call error_mesg ('sealw99_mod', &
'initial '//trim(gas)//' transmission function&
& calculation uses '//trim(gas)//' mixing &
&ratio of:' // chvers7 //', which is the value &
&at the current time.', NOTE)
endif
endif
endif
endif !(alarm == 0)
!---------------------------------------------------------------------
! if it is desired that the tfs be calculated only on the first
! step, call the appropriate subroutines to do so. redefine the
! gas values used for the last tf calculation redefine the
! gas values used for the last tf calculation.
!---------------------------------------------------------------------
! else ! (.not. calc_gas_tfs_on_first_step)
else if (calc_gas_tfs_on_first_step) then
if (trim(gas) == 'ch4') then
call ch4_time_vary (gas_for_next_tf_calc)
write (chvers7, '(4pe15.7)') gas_for_next_tf_calc
else if (trim(gas) == 'n2o') then
call n2o_time_vary (gas_for_next_tf_calc)
write (chvers7, '(3pe15.7)') gas_for_next_tf_calc
else if (trim(gas) == 'co2') then
call co2_time_vary (gas_for_next_tf_calc)
write (chvers7, '(3pe15.7)') gas_for_next_tf_calc
endif
gas_for_last_tf_calc = gas_for_next_tf_calc
!---------------------------------------------------------------------
! if a record of the tf calculation path is desired, print out the
! relevant data.
!---------------------------------------------------------------------
if (verbose >= 1) then
if (gas_tf_offset /= 0.0) then
call error_mesg ('sealw99_mod', &
'calculating '//trim(gas)//' transmission functions&
& at time ' //chvers//chvers2//chvers3//chvers4// &
chvers5//chvers6// ', using '//trim(gas)//' mixing &
&ratio of:' // chvers7 //', which is the value ' &
//chvers10// 'hours from current time.', NOTE)
else
call error_mesg ('sealw99_mod', &
'calculating '//trim(gas)//' transmission functions&
& at time ' //chvers//chvers2//chvers3//chvers4//&
chvers5//chvers6// ', using '//trim(gas)//' mixing &
&ratio of:' // chvers7 //', which is the value at &
&the current time.', NOTE)
endif
endif
!---------------------------------------------------------------------
! set the flag to indicate that the initial tf calculation has been
! completed.
!---------------------------------------------------------------------
do_gas_tf_calc_init = .false.
else if (calc_gas_tfs_monthly) then
if (do_gas_tf_calc_init) then
if (trim(gas) == 'ch4') then
call ch4_time_vary (gas_for_next_tf_calc)
write (chvers7, '(4pe15.7)') gas_for_next_tf_calc
else if (trim(gas) == 'n2o') then
call n2o_time_vary (gas_for_next_tf_calc)
write (chvers7, '(3pe15.7)') gas_for_next_tf_calc
else if (trim(gas) == 'co2') then
call co2_time_vary (gas_for_next_tf_calc)
write (chvers7, '(3pe15.7)') gas_for_next_tf_calc
endif
gas_for_last_tf_calc = gas_for_next_tf_calc
!---------------------------------------------------------------------
! if a record of the tf calculation path is desired, print out the
! relevant data.
!---------------------------------------------------------------------
if (verbose >= 1) then
if (gas_tf_offset /= 0.0) then
call error_mesg ('sealw99_mod', &
'calculating '//trim(gas)//' transmission functions&
& at time ' //chvers//chvers2//chvers3//chvers4// &
chvers5//chvers6// ', using '//trim(gas)//' mixing &
&ratio of:' // chvers7 //', which is the value ' &
//chvers10// 'hours from current time .', NOTE)
else
call error_mesg ('sealw99_mod', &
'calculating '//trim(gas)//' transmission functions&
& at time ' //chvers//chvers2//chvers3//chvers4//&
chvers5//chvers6// ', using '//trim(gas)//' mixing &
&ratio of:' // chvers7 //', which is the value at &
&the current time.', NOTE)
endif
endif
!---------------------------------------------------------------------
! set the flag to indicate that the initial tf calculation has been
! completed.
!---------------------------------------------------------------------
do_gas_tf_calc_init = .false.
month_of_gas_tf_calc = month
else ! (do_gas_tf_calc_init)
if (month /= month_of_gas_tf_calc) then
if (trim(gas) == 'ch4') then
call ch4_time_vary (gas_for_next_tf_calc)
write (chvers7, '(4pe15.7)') gas_for_next_tf_calc
else if (trim(gas) == 'n2o') then
call n2o_time_vary (gas_for_next_tf_calc)
write (chvers7, '(3pe15.7)') gas_for_next_tf_calc
else if (trim(gas) == 'co2') then
call co2_time_vary (gas_for_next_tf_calc)
write (chvers7, '(3pe15.7)') gas_for_next_tf_calc
endif
gas_for_last_tf_calc = gas_for_next_tf_calc
!---------------------------------------------------------------------
! if a record of the tf calculation path is desired, print out the
! relevant data.
!---------------------------------------------------------------------
if (verbose >= 1) then
if (gas_tf_offset /= 0.0) then
call error_mesg ('sealw99_mod', &
'calculating '//trim(gas)//' transmission functions&
& at time ' //chvers//chvers2//chvers3//chvers4// &
chvers5//chvers6// ', using '//trim(gas)//' mixing &
&ratio of:' // chvers7 //', which is the value ' &
//chvers10// 'hours from current time.', NOTE)
else
call error_mesg ('sealw99_mod', &
'calculating '//trim(gas)//' transmission functions&
& at time ' //chvers//chvers2//chvers3//chvers4//&
chvers5//chvers6// ', using '//trim(gas)//' mixing &
&ratio of:' // chvers7 //', which is the value at &
&the current time.', NOTE)
endif
endif
!---------------------------------------------------------------------
! set the flag to indicate that the initial tf calculation has been
! completed.
!---------------------------------------------------------------------
month_of_gas_tf_calc = month
endif
endif ! (do_gas_tf_calc_init)
endif ! (.not. calc_gas_tfs_on_first_step)
!---------------------------------------------------------------------
! set the flag to indicate that the tf calculation has been
! completed on the current timestep.
!---------------------------------------------------------------------
do_gas_tf_calc = .false.
!---------------------------------------------------------------------
! if the time variation of the gas has not yet begun, and it is the
! initial call of the job, call the appropriate subroutines to
! define the transmission functions.
!---------------------------------------------------------------------
else ! (Rad_time >= Gas_time)
if (do_gas_tf_calc_init) then
if (trim(gas) == 'ch4') then
call ch4_time_vary (gas_for_last_tf_calc)
write (chvers7, '(4pe15.7)') gas_for_last_tf_calc
else if (trim(gas) == 'n2o') then
call n2o_time_vary (gas_for_last_tf_calc)
write (chvers7, '(3pe15.7)') gas_for_last_tf_calc
else if (trim(gas) == 'co2') then
call co2_time_vary (gas_for_last_tf_calc)
write (chvers7, '(3pe15.7)') gas_for_last_tf_calc
endif
!---------------------------------------------------------------------
! set the flag to indicate that the initial tf calculation has been
! completed.
!---------------------------------------------------------------------
do_gas_tf_calc_init = .false.
!---------------------------------------------------------------------
! if a record of the tf calculation path is desired, print out the
! relevant data.
!---------------------------------------------------------------------
if (verbose >= 1) then
write (chvers10, '( f9.3)') gas_tf_offset
if (gas_tf_offset /= 0.0) then
call error_mesg ('sealw99_mod', &
'initial '//trim(gas)//' transmission function &
&calculation uses '//trim(gas)//' mixing ratio of:' &
// chvers7 //', which is the value '//chvers10// &
'hours from current time.', NOTE)
else
call error_mesg ('sealw99_mod', &
'initial '//trim(gas)//' transmission function &
&calculation uses '//trim(gas)//' mixing ratio of:' &
// chvers7 //', which is the value at the current &
&time.', NOTE)
endif
endif
endif ! (do_gas_tf_calc_init)
endif ! (Rad_time >= Gas_time)
!--------------------------------------------------------------------
end subroutine obtain_gas_tfs
!#####################################################################
!
!
! Subroutine to allocate variables needed for longwave diagnostics
!
!
! call sealw99_alloc (ix, jx, kx, Lw_diagnostics)
!
!
! Dimension 1 length of radiation arrays to be allocated
!
!
! Dimension 2 length of radiation arrays to be allocated
!
!
! Dimension 3 length of radiation arrays to be allocated
!
!
! lw_diagnostics_type variable containing longwave
! radiation output data
!
!
!
subroutine sealw99_alloc (ix, jx, kx, Lw_diagnostics)
!--------------------------------------------------------------------
! sealw99_alloc allocates and initializes the components of the
! lw_diagnostics_type variable Lw_diagnostics which holds diagnostic
! output generated by sealw99_mod.
!--------------------------------------------------------------------
integer, intent(in) :: ix, jx, kx
type(lw_diagnostics_type), intent(inout) :: Lw_diagnostics
!--------------------------------------------------------------------
! intent(in) variables:
!
! ix,jx,kx (i,j,k) lengths of radiation arrays to be allocated
!
!
! intent(inout) variables:
!
! Lw_diagnostics
! lw_diagnostics_type variable containing diagnostic
! longwave output used by the radiation diagnostics
! module
!
!---------------------------------------------------------------------
!--------------------------------------------------------------------
! local variables:
integer :: NBTRGE, NBLY
!---------------------------------------------------------------------
! allocate (and initialize where necessary) lw_diagnostics_type
! component arrays.
!---------------------------------------------------------------------
NBTRGE = Lw_parameters%NBTRGE
NBLY = Lw_parameters%NBLY
allocate ( Lw_diagnostics%flx1e1 (ix, jx ) )
allocate ( Lw_diagnostics%fluxn (ix, jx, kx+1, 6+NBTRGE) )
allocate (Lw_diagnostics%cts_out (ix, jx, kx, 6 ) )
allocate (Lw_diagnostics%cts_outcf (ix, jx, kx, 6 ) )
allocate (Lw_diagnostics%gxcts (ix, jx ) )
allocate (Lw_diagnostics%excts (ix, jx, kx ) )
allocate (Lw_diagnostics%exctsn (ix, jx, kx, NBLY ) )
allocate (Lw_diagnostics%fctsg (ix, jx, NBLY ) )
Lw_diagnostics%flx1e1 = 0.
Lw_diagnostics%cts_out = 0.
Lw_diagnostics%cts_outcf = 0.
Lw_diagnostics%gxcts = 0.
Lw_diagnostics%excts = 0.
Lw_diagnostics%exctsn = 0.
Lw_diagnostics%fctsg = 0.
Lw_diagnostics%fluxn (:,:,:,:) = 0.0
if (Rad_control%do_totcld_forcing) then
allocate ( Lw_diagnostics%fluxncf (ix, jx, kx+1, 6+NBTRGE) )
Lw_diagnostics%fluxncf(:,:,:,:) = 0.0
endif
allocate( Lw_diagnostics%flx1e1f (ix, jx, NBTRGE ) )
Lw_diagnostics%flx1e1f = 0.
!--------------------------------------------------------------------
end subroutine sealw99_alloc
!####################################################################
!
!
! Subroutine the calculate the cool to space approximation longwave
! radiation.
!
!
! call cool_to_space_approx ( pflux_in, source, &
! trans, cld_trans, cld_ind, &
! Lw_diagnostics, &
! trans2 )
!
!
! pressure values at flux levels
!
!
! band integrated longwave source function of each model layer
!
!
! clear sky longwave transmission
!
!
! cloud transmission
!
!
! cloud type index
!
!
! longwave diagnostics output
!
!
! optional input alternative transmission profile
!
!
!
subroutine cool_to_space_approx ( pflux_in, source, trans, cld_trans, &
cld_ind, Lw_diagnostics, trans2 )
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
real, dimension (:,:,:), intent(in) :: pflux_in
real, dimension (:,:,:,:), intent(in) :: source, trans, &
cld_trans
integer, dimension (:), intent(in) :: cld_ind
type(lw_diagnostics_type), intent(inout) :: Lw_diagnostics
real, dimension (:,:,:), intent(in), optional :: trans2
!---------------------------------------------------------------------
! intent(in) variables:
!
! pflux_in
! source
! trans
! cld_trans
! cld_ind
!
! intent(inout) variables:
!
! Lw_diagnostics
!
! intent(in),optional:
! trans2
!
!---------------------------------------------------------------------
!---------------------------------------------------------------------
! local variables:
real, dimension(size(pflux_in,1), size(pflux_in,2), &
size(pflux_in,3)-1) :: pdfinv
integer :: i,j,kk
integer :: index, nbands
!---------------------------------------------------------------------
! local variables:
!
! pdfinv
! index
! nbands
! j
!
!---------------------------------------------------------------------
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
nbands = size(source,4) - NBTRGE
do kk = KS,KE
do j = 1,size(pdfinv(:,:,:),2)
do i = 1,size(pdfinv(:,:,:),1)
pdfinv(i,j,kk) = 1.0/(pflux_in(i,j,kk+KS+1-(KS)) - pflux_in(i,j,kk))
end do
end do
end do
!---------------------------------------------------------------------
!
!--------------------------------------------------------------------
do index=1,nbands
if (index == 1 ) then
do kk = KS,KE
do j = 1,size(Lw_diagnostics%cts_out(:,:,:,:),2)
do i = 1,size(Lw_diagnostics%cts_out(:,:,:,:),1)
Lw_diagnostics%cts_out(i,j,kk,index) = (radcon_mks*pdfinv(i,j,kk)* &
source(i,j,kk,index)* &
(trans(i,j,kk,index)* &
cld_trans(i,j,kk+KS+1-(KS), cld_ind(index)) - &
trans2(i,j,kk)* &
cld_trans(i,j,kk, cld_ind(index)))*1.0e-03)
end do
end do
end do
else
do kk = KS,KE
do j = 1,size(Lw_diagnostics%cts_out(:,:,:,:),2)
do i = 1,size(Lw_diagnostics%cts_out(:,:,:,:),1)
Lw_diagnostics%cts_out(i,j,kk,index) = (radcon_mks*pdfinv(i,j,kk)* &
source(i,j,kk, index)* &
(trans(i,j,kk+KS+1-(KS),index )* &
cld_trans(i,j,kk+KS+1-(KS), cld_ind(index)) - &
trans(i,j,kk,index )* &
cld_trans(i,j,kk, cld_ind(index)))*1.0e-03)
end do
end do
end do
endif
if (Rad_control%do_totcld_forcing) then
if (index == 1 ) then
do kk = KS,KE
do j = 1,size(Lw_diagnostics%cts_outcf(:,:,:,:),2)
do i = 1,size(Lw_diagnostics%cts_outcf(:,:,:,:),1)
Lw_diagnostics%cts_outcf(i,j,kk,index) = (radcon_mks*pdfinv(i,j,kk)* &
source(i,j,kk,index)* &
(trans(i,j,kk,index ) - &
trans2(i,j,kk))*1.0e-03)
end do
end do
end do
else
do kk = KS,KE
do j = 1,size(Lw_diagnostics%cts_outcf(:,:,:,:),2)
do i = 1,size(Lw_diagnostics%cts_outcf(:,:,:,:),1)
Lw_diagnostics%cts_outcf(i,j,kk,index) = (radcon_mks*pdfinv(i,j,kk)* &
source(i,j,kk,index)* &
(trans(i,j,kk+KS+1-(KS), index ) - &
trans(i,j,kk, index ))*1.0e-03)
end do
end do
end do
endif
endif
end do ! (index loop)
!--------------------------------------------------------------------
end subroutine cool_to_space_approx
!####################################################################
!
!
! cool_to_space calculates the cool-to-space cooling rate for
! a band n.
!
!
! call cool_to_space_exact ( cldtf, &
! Atmos_input, Optical, Gas_tf, &
! sorc, to3cnt, Lw_diagnostics, &
! cts_sum, cts_sumcf, &
! gxctscf)
!
!
! cloud transmission function between levels k level KS.
!
!
! Atmospheric input to the cool to space approximation method
!
!
! Optical depth of atmospheric layers and clouds
!
!
! Gas transmission function
!
!
! band-integrated Planck function, for each combined
! band in the 160-1200 cm-1 region.
!
!
! transmission functions between levels k and
! level KS for the 990-1070 cm-1 range.
!
!
! Longwave diagnostics
!
!
! Cool to space heating rates
!
!
! Cool to space heating rates due to cloud forcing
!
!
! gxcts is the "exact" surface flux accumulated over
! the frequency bands in the 160-1200 cm-1 range.
!
!
subroutine cool_to_space_exact (cldtf, Atmos_input, Optical, Gas_tf, &
sorc, to3cnt, Lw_diagnostics, &
cts_sum, cts_sumcf, gxctscf, &
including_aerosols)
!-----------------------------------------------------------------------
! cool_to_space calculates the cool-to-space cooling rate for
! a band n.
!-----------------------------------------------------------------------
real, dimension (:,:,:,:), intent(in) :: cldtf
type(atmos_input_type), intent(in) :: Atmos_input
type(optical_path_type), intent(inout) :: Optical
type(gas_tf_type), intent(inout) :: Gas_tf
real, dimension (:,:,:,:), intent(in) :: sorc
real, dimension (:,:,:), intent(in) :: to3cnt
type(lw_diagnostics_type), intent(inout) :: Lw_diagnostics
real, dimension(:,:,:), intent(inout) :: cts_sum, cts_sumcf
real, dimension(:,:), intent(inout) :: gxctscf
logical, intent(in) :: including_aerosols
!--------------------------------------------------------------------
! intent(in) variables:
!
! cldtf cloud transmission function between levels k and
! level KS.
! Atmos_input
! sorc band-integrated Planck function, for each combined
! band in the 160-1200 cm-1 region.
! to3cnt transmission functions between levels k and
! level KS for the 990-1070 cm-1 range.
!
! intent(inout) variables:
!
! Optical
! Gas_tf
! Lw_diagnostics
! cts_sum
! cts_sumcf
! gxctscf
!
!--------------------------------------------------------------------
!-----------------------------------------------------------------------
! local variables
!-----------------------------------------------------------------------
integer :: i,kk
real, dimension(size(Atmos_input%pflux,1), &
size(Atmos_input%pflux,2), &
size(Atmos_input%pflux,3)-1) :: pdfinv, pdfinv2
real, dimension(size(Atmos_input%pflux,1), &
size(Atmos_input%pflux,2), &
size(Atmos_input%pflux,3) ) :: &
dte1, press, temp, pflux
integer, dimension(size(Atmos_input%pflux,1), &
size(Atmos_input%pflux,2), &
size(Atmos_input%pflux,3) ) :: ixoe1
real, dimension(size(Atmos_input%pflux,1), &
size(Atmos_input%pflux,2)) :: &
pfac1, pfac2
real, dimension(size(Atmos_input%pflux,1), &
size(Atmos_input%pflux,2), &
size(Atmos_input%pflux,3)) :: &
sorc_tmp, ctmp, totch2o_tmp, totaer_tmp
real, dimension(size(Atmos_input%pflux,1), &
size(Atmos_input%pflux,2), &
2:size(Atmos_input%pflux,3)) :: &
totvo2_tmp
real, dimension(size(Atmos_input%pflux,1), &
size(Atmos_input%pflux,2), &
size(Atmos_input%pflux,3)-1) :: &
exctscf, tt, x, y, topm, &
topphi, phitmp, psitmp, ag, &
agg, f, ff, tmp1, tmp2, fac1, &
fac2, cfc_tf
real, dimension(size(Atmos_input%pflux,1), &
size(Atmos_input%pflux,2), &
size(Atmos_input%pflux,3)-1, NBLY) :: &
exctsncf
real, dimension(size(Atmos_input%pflux,1), &
size(Atmos_input%pflux,2), &
NBLY) :: &
fctsgcf
integer :: n, k, j, ioffset
!-----------------------------------------------------------------------
! local variables
!
! pdfinv inverse of pressure difference between flux levels.
! pdfinv2
! dte1
! press pressure at data levels of model.
! temp temperature at data levels of model.
! pflux pressure at flux levels of model.
! ixoe1
! pfac1
! pfac2
! sorc_tmp
! ctmp
! totch2o_tmp
! totaer_tmp
! totvo2_tmp
! exctscf
! tt
! x
! y
! topm
! topphi
! phitmp
! psitmp
! ag
! agg
! f
! ff
! tmp1
! tmp2
! fac1
! fac2
! cfc_tf
!
!---------------------------------------------------------------------
!---------------------------------------------------------------------
! convert press and pflux to cgs.
do kk = 1,size(press(:,:,:),3)
do j = 1,size(press(:,:,:),2)
do i = 1,size(press(:,:,:),1)
press(i,j,kk) = 10.0*Atmos_input%press(i,j,kk)
end do
end do
end do
do kk = 1,size(pflux(:,:,:),3)
do j = 1,size(pflux(:,:,:),2)
do i = 1,size(pflux(:,:,:),1)
pflux(i,j,kk) = 10.0*Atmos_input%pflux(i,j,kk)
end do
end do
end do
do kk = 1,size(temp(:,:,:),3)
do j = 1,size(temp(:,:,:),2)
do i = 1,size(temp(:,:,:),1)
temp(i,j,kk) = Atmos_input%temp(i,j,kk)
end do
end do
end do
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
do kk = KS,KE
do j = 1,size(pdfinv2(:,:,:),2)
do i = 1,size(pdfinv2(:,:,:),1)
pdfinv2(i,j,kk) = 1.0/(pflux(i,j,kk+KS+1-(KS)) - pflux(i,j,kk))
end do
end do
end do
do kk = KS,KE
do j = 1,size(pdfinv(:,:,:),2)
do i = 1,size(pdfinv(:,:,:),1)
pdfinv(i,j,kk) = 1.0/(Atmos_input%pflux(i,j,kk+KS+1-(KS)) - &
Atmos_input%pflux(i,j,kk))
end do
end do
end do
!----------------------------------------------------------------------
ioffset = Lw_parameters%offset
!-----------------------------------------------------------------------
! initialize quantities.
!-----------------------------------------------------------------------
do kk = KS,KE
do j = 1,size(Lw_diagnostics%excts(:,:,:),2)
do i = 1,size(Lw_diagnostics%excts(:,:,:),1)
Lw_diagnostics%excts(i,j,kk) = 0.0E+00
end do
end do
end do
do j = 1,size(Lw_diagnostics%gxcts(:,:),2)
do i = 1,size(Lw_diagnostics%gxcts(:,:),1)
Lw_diagnostics%gxcts(i,j) = 0.0E+00
end do
end do
if (Rad_control%do_totcld_forcing) then
do kk = KS,KE
do j = 1,size(exctscf(:,:,:),2)
do i = 1,size(exctscf(:,:,:),1)
exctscf(i,j,kk) = 0.0E+00
end do
end do
end do
do j = 1,size(gxctscf(:,:),2)
do i = 1,size(gxctscf(:,:),1)
gxctscf(i,j) = 0.0E+00
end do
end do
endif
!-----------------------------------------------------------------------
! compute temperature quantities.
!-----------------------------------------------------------------------
do kk = KS,KE
do j = 1,size(x(:,:,:),2)
do i = 1,size(x(:,:,:),1)
x(i,j,kk) = temp(i,j,kk) - 2.5E+02
end do
end do
end do
do kk = KS,KE
do j = 1,size(y(:,:,:),2)
do i = 1,size(y(:,:,:),1)
y(i,j,kk) = x(i,j,kk)*x(i,j,kk)
end do
end do
end do
do j = 1,size(ctmp(:,:,:),2)
do i = 1,size(ctmp(:,:,:),1)
ctmp(i,j,KS) = 1.0E+00
end do
end do
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
call locate_in_table(temp_1, temp, dte1, ixoe1, KS, KE+1)
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
do j = 1,size(Lw_diagnostics%fctsg(:,:,:),2)
do i = 1,size(Lw_diagnostics%fctsg(:,:,:),1)
Lw_diagnostics%fctsg(i,j,NBLY) = 0.0
end do
end do
do n=1,NBLY-1
!-----------------------------------------------------------------------
! obtain temperature correction capphi, cappsi, then multiply
! by optical path var1, var2 to compute temperature-corrected
! optical path and mean pressure for a layer: phitmp, psitmp.
!-----------------------------------------------------------------------
do kk = KS,KE
do j = 1,size(f(:,:,:),2)
do i = 1,size(f(:,:,:),1)
f(i,j,kk) = 0.44194E-01*(apcm (n)*x(i,j,kk) + &
bpcm (n)*y(i,j,kk))
end do
end do
end do
do kk = KS,KE
do j = 1,size(ff(:,:,:),2)
do i = 1,size(ff(:,:,:),1)
ff(i,j,kk) = 0.44194E-01*(atpcm(n)*x(i,j,kk) + &
btpcm(n)*y(i,j,kk))
end do
end do
end do
do kk = KS,KE
do j = 1,size(ag(:,:,:),2)
do i = 1,size(ag(:,:,:),1)
ag(i,j,kk) = (1.418191E+00 + f (i,j,kk))* &
f (i,j,kk) + 1.0E+00
end do
end do
end do
do kk = KS,KE
do j = 1,size(agg(:,:,:),2)
do i = 1,size(agg(:,:,:),1)
agg(i,j,kk) = (1.418191E+00 + ff(i,j,kk))* &
ff(i,j,kk) + 1.0E+00
end do
end do
end do
do kk = KS,KE
do j = 1,size(phitmp(:,:,:),2)
do i = 1,size(phitmp(:,:,:),1)
phitmp(i,j,kk) = Optical%var1(i,j,kk)* &
((((ag (i,j,kk)* &
ag (i,j,kk))**2)**2)**2)
end do
end do
end do
do kk = KS,KE
do j = 1,size(psitmp(:,:,:),2)
do i = 1,size(psitmp(:,:,:),1)
psitmp(i,j,kk) = Optical%var2(i,j,kk)* &
((((agg(i,j,kk)* &
agg(i,j,kk))**2)**2)**2)
end do
end do
end do
!-----------------------------------------------------------------------
! obtain optical path and mean pressure from the top of the
! atmosphere to the level k.
!-----------------------------------------------------------------------
do j = 1,size(topm(:,:,:),2)
do i = 1,size(topm(:,:,:),1)
topm(i,j,KS) = phitmp(i,j,KS)
end do
end do
do j = 1,size(topphi(:,:,:),2)
do i = 1,size(topphi(:,:,:),1)
topphi(i,j,KS) = psitmp(i,j,KS)
end do
end do
do k=KS+1,KE
do j = 1,size(topm(:,:,:),2)
do i = 1,size(topm(:,:,:),1)
topm(i,j,k) = topm (i,j,k-1) + phitmp(i,j,k)
end do
end do
do j = 1,size(topphi(:,:,:),2)
do i = 1,size(topphi(:,:,:),1)
topphi(i,j,k) = topphi(i,j,k-1) + psitmp(i,j,k)
end do
end do
enddo
!-----------------------------------------------------------------------
! tt is the cloud-free h2o cool-to-space transmission function.
!-----------------------------------------------------------------------
if (Lw_control%do_h2o) then
do kk = KS,KE
do j = 1,size(fac1(:,:,:),2)
do i = 1,size(fac1(:,:,:),1)
fac1(i,j,kk) = acomb(n)*topm(i,j,kk)
end do
end do
end do
do kk = KS,KE
do j = 1,size(fac2(:,:,:),2)
do i = 1,size(fac2(:,:,:),1)
fac2(i,j,kk) = fac1(i,j,kk)*topm(i,j,kk)/ &
(bcomb(n)*topphi(i,j,kk))
end do
end do
end do
do kk = KS,KE
do j = 1,size(tmp1(:,:,:),2)
do i = 1,size(tmp1(:,:,:),1)
tmp1(i,j,kk) = fac1(i,j,kk)/SQRT(1.0E+00 + &
fac2(i,j,kk))
end do
end do
end do
else
do kk = KS,KE
do j = 1,size(tmp1(:,:,:),2)
do i = 1,size(tmp1(:,:,:),1)
tmp1(i,j,kk) = 0.0
end do
end do
end do
endif
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
if (n >= band_no_start(1) .and. n <= band_no_end(1)) then
! 160-400 cm-1 region (h2o)
if (trim(Lw_control%continuum_form) == 'ckd2.1' .or. &
trim(Lw_control%continuum_form) == 'ckd2.4' ) then
! bands 1-24.
call get_totch2o (n, Optical, totch2o_tmp, dte1, ixoe1)
do kk = KS,KE
do j = 1,size(tt(:,:,:),2)
do i = 1,size(tt(:,:,:),1)
tt(i,j,kk) = EXP(-1.0*(tmp1(i,j,kk) + diffac* &
totch2o_tmp(i,j,kk+KS+1-(KS))))
end do
end do
end do
else if (trim(Lw_control%continuum_form) == 'rsb' ) then
! bands 1-4.
do kk = KS,KE
do j = 1,size(tt(:,:,:),2)
do i = 1,size(tt(:,:,:),1)
tt(i,j,kk) = EXP(-1.0*tmp1(i,j,kk))
end do
end do
end do
endif
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
else if (n >= band_no_start(2) .and. n <= band_no_end(2)) then
! 400-560 cm-1 region (h2o)
if (trim(Lw_control%continuum_form) == 'ckd2.1' .or. &
trim(Lw_control%continuum_form) == 'ckd2.4' ) then
! bands 25-40.
call get_totch2o (n, Optical, totch2o_tmp, dte1, ixoe1)
do kk = KS,KE
do j = 1,size(tt(:,:,:),2)
do i = 1,size(tt(:,:,:),1)
tt(i,j,kk) = EXP(-1.0*(tmp1(i,j,kk) + diffac* &
totch2o_tmp(i,j,kk+KS+1-(KS))))
end do
end do
end do
else if (trim(Lw_control%continuum_form) == 'rsb' ) then
! bands 5-8.
call get_totvo2 (n, Optical, totvo2_tmp)
do kk = KS,KE
do j = 1,size(tt(:,:,:),2)
do i = 1,size(tt(:,:,:),1)
tt(i,j,kk) = EXP(-1.0*(tmp1(i,j,kk) + &
totvo2_tmp(i,j,kk+KS+1-(KS))))
end do
end do
end do
endif
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
else if (n >= band_no_start(3) .and. n <= band_no_end(3)) then
! 560-630 cm-1 region (h2o, co2, n2o, aerosol)
if (trim(Lw_control%continuum_form) == 'ckd2.1' .or. &
trim(Lw_control%continuum_form) == 'ckd2.4' ) then
! band 41.
call get_totch2obd (n-40, Optical, totch2o_tmp)
do kk = KS,KE
do j = 1,size(tmp2(:,:,:),2)
do i = 1,size(tmp2(:,:,:),1)
tmp2(i,j,kk) = tmp1(i,j,kk) + diffac* &
totch2o_tmp(i,j,kk+KS+1-(KS))
end do
end do
end do
else if (trim(Lw_control%continuum_form) == 'rsb' ) then
! band 9
call get_totvo2 (n, Optical, totvo2_tmp)
do kk = KS,KE
do j = 1,size(tmp2(:,:,:),2)
do i = 1,size(tmp2(:,:,:),1)
tmp2(i,j,kk) = tmp1(i,j,kk) + &
totvo2_tmp(i,j,kk+KS+1-(KS))
end do
end do
end do
endif
if (including_aerosols) then
do kk = 1,size(totaer_tmp(:,:,:),3)
do j = 1,size(totaer_tmp(:,:,:),2)
do i = 1,size(totaer_tmp(:,:,:),1)
totaer_tmp(i,j,kk) = Optical%totaerooptdep(i,j,kk,1)
end do
end do
end do
do kk = KS,KE
do j = 1,size(tmp2(:,:,:),2)
do i = 1,size(tmp2(:,:,:),1)
tmp2(i,j,kk) = tmp2(i,j,kk) + &
totaer_tmp (i,j,kk+KS+1-(KS))
end do
end do
end do
endif
do kk = KS,KE
do j = 1,size(tt(:,:,:),2)
do i = 1,size(tt(:,:,:),1)
tt(i,j,kk) = EXP(-1.0E+00*tmp2(i,j,kk))* &
(Gas_tf%co2spnb(i,j,kk+KS+1-(KS),1)* &
Gas_tf%tn2o17(i,j,kk+KS+1-(KS)))
end do
end do
end do
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
else if (n >= band_no_start(4) .and. n <= band_no_end(4)) then
! 630-700 cm-1 region (h2o, co2, aerosol)
if (trim(Lw_control%continuum_form) == 'ckd2.1' .or. &
trim(Lw_control%continuum_form) == 'ckd2.4' ) then
! band 42.
call get_totch2obd (n-40, Optical, totch2o_tmp)
do kk = KS,KE
do j = 1,size(tmp2(:,:,:),2)
do i = 1,size(tmp2(:,:,:),1)
tmp2(i,j,kk) = tmp1(i,j,kk) + diffac* &
totch2o_tmp(i,j,kk+KS+1-(KS))
end do
end do
end do
else if (trim(Lw_control%continuum_form) == 'rsb' ) then
! band 10
call get_totvo2 (n, Optical, totvo2_tmp)
do kk = KS,KE
do j = 1,size(tmp2(:,:,:),2)
do i = 1,size(tmp2(:,:,:),1)
tmp2(i,j,kk) = tmp1(i,j,kk) + &
totvo2_tmp(i,j,kk+KS+1-(KS))
end do
end do
end do
endif
if (including_aerosols) then
do kk = 1,size(totaer_tmp(:,:,:),3)
do j = 1,size(totaer_tmp(:,:,:),2)
do i = 1,size(totaer_tmp(:,:,:),1)
totaer_tmp(i,j,kk) = Optical%totaerooptdep(i,j,kk,2)
end do
end do
end do
do kk = KS,KE
do j = 1,size(tmp2(:,:,:),2)
do i = 1,size(tmp2(:,:,:),1)
tmp2(i,j,kk) = tmp2(i,j,kk) + &
totaer_tmp (i,j,kk+KS+1-(KS))
end do
end do
end do
endif
do kk = KS,KE
do j = 1,size(tt(:,:,:),2)
do i = 1,size(tt(:,:,:),1)
tt(i,j,kk) = EXP(-1.0E+00*tmp2(i,j,kk))* &
Gas_tf%co2spnb(i,j,kk+KS+1-(KS),2)
end do
end do
end do
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
else if (n >= band_no_start(5) .and. n <= band_no_end(5)) then
! 700-800 cm-1 region (h2o, co2, aerosol)
if (trim(Lw_control%continuum_form) == 'ckd2.1' .or. &
trim(Lw_control%continuum_form) == 'ckd2.4' ) then
! band 43.
call get_totch2obd (n-40, Optical, totch2o_tmp)
do kk = KS,KE
do j = 1,size(tmp2(:,:,:),2)
do i = 1,size(tmp2(:,:,:),1)
tmp2(i,j,kk) = tmp1(i,j,kk) + diffac* &
totch2o_tmp(i,j,kk+KS+1-(KS))
end do
end do
end do
else if (trim(Lw_control%continuum_form) == 'rsb' ) then
! band 11
call get_totvo2 (n, Optical, totvo2_tmp)
do kk = KS,KE
do j = 1,size(tmp2(:,:,:),2)
do i = 1,size(tmp2(:,:,:),1)
tmp2(i,j,kk) = tmp1(i,j,kk) + &
totvo2_tmp(i,j,kk+KS+1-(KS))
end do
end do
end do
endif
if (including_aerosols) then
do kk = 1,size(totaer_tmp(:,:,:),3)
do j = 1,size(totaer_tmp(:,:,:),2)
do i = 1,size(totaer_tmp(:,:,:),1)
totaer_tmp(i,j,kk) = Optical%totaerooptdep(i,j,kk,3)
end do
end do
end do
do kk = KS,KE
do j = 1,size(tmp2(:,:,:),2)
do i = 1,size(tmp2(:,:,:),1)
tmp2(i,j,kk) = tmp2(i,j,kk) + &
totaer_tmp (i,j,kk+KS+1-(KS))
end do
end do
end do
endif
do kk = KS,KE
do j = 1,size(tt(:,:,:),2)
do i = 1,size(tt(:,:,:),1)
tt(i,j,kk) = EXP(-1.0E+00*tmp2(i,j,kk))* &
Gas_tf%co2spnb(i,j,kk+KS+1-(KS),3)
end do
end do
end do
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
else if (n >= band_no_start(6) .and. n <= band_no_end(6)) then
! 800-990 cm-1 region (h2o, aerosol)
if (trim(Lw_control%continuum_form) == 'ckd2.1' .or. &
trim(Lw_control%continuum_form) == 'ckd2.4' ) then
! bands 44-45.
call get_totch2obd (n-40, Optical, totch2o_tmp)
do kk = KS,KE
do j = 1,size(tmp2(:,:,:),2)
do i = 1,size(tmp2(:,:,:),1)
tmp2(i,j,kk) = tmp1(i,j,kk) + diffac* &
totch2o_tmp(i,j,kk+KS+1-(KS))
end do
end do
end do
else if (trim(Lw_control%continuum_form) == 'rsb' ) then
! bands 12-13
call get_totvo2 (n, Optical, totvo2_tmp)
do kk = KS,KE
do j = 1,size(tmp2(:,:,:),2)
do i = 1,size(tmp2(:,:,:),1)
tmp2(i,j,kk) = tmp1(i,j,kk) + &
totvo2_tmp(i,j,kk+KS+1-(KS))
end do
end do
end do
endif
if (including_aerosols) then
do kk = 1,size(totaer_tmp(:,:,:),3)
do j = 1,size(totaer_tmp(:,:,:),2)
do i = 1,size(totaer_tmp(:,:,:),1)
totaer_tmp(i,j,kk) = Optical%totaerooptdep(i,j,kk,n-8-ioffset)
end do
end do
end do
do kk = KS,KE
do j = 1,size(tmp2(:,:,:),2)
do i = 1,size(tmp2(:,:,:),1)
tmp2(i,j,kk) = tmp2(i,j,kk) + &
totaer_tmp (i,j,kk+KS+1-(KS))
end do
end do
end do
endif
do kk = KS,KE
do j = 1,size(tt(:,:,:),2)
do i = 1,size(tt(:,:,:),1)
tt(i,j,kk) = EXP(-1.0E+00*tmp2(i,j,kk))
end do
end do
end do
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
else if (n >= band_no_start(7) .and. n <= band_no_end(7)) then
! 990-1070 cm-1 region (h2o(lines), o3)
! band 46 (ckd2.1) or 14 (rsb)
do kk = KS,KE
do j = 1,size(tt(:,:,:),2)
do i = 1,size(tt(:,:,:),1)
tt(i,j,kk) = EXP(-1.0E+00*tmp1(i,j,kk))* &
to3cnt (i,j,kk+KS+1-(KS))
end do
end do
end do
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
else if (n >= band_no_start(8) .and. n <= band_no_end(8)) then
! 1070-1200 cm-1 region (h2o, n2o)
if (trim(Lw_control%continuum_form) == 'ckd2.1' .or. &
trim(Lw_control%continuum_form) == 'ckd2.4' ) then
! band 47.
call get_totch2obd (n-40, Optical, totch2o_tmp)
do kk = KS,KE
do j = 1,size(tmp2(:,:,:),2)
do i = 1,size(tmp2(:,:,:),1)
tmp2(i,j,kk) = tmp1(i,j,kk) + diffac* &
totch2o_tmp(i,j,kk+KS+1-(KS))
end do
end do
end do
else if (trim(Lw_control%continuum_form) == 'rsb' ) then
! band 15
call get_totvo2 (n, Optical, totvo2_tmp)
do kk = KS,KE
do j = 1,size(tmp2(:,:,:),2)
do i = 1,size(tmp2(:,:,:),1)
tmp2(i,j,kk) = tmp1(i,j,kk) + &
totvo2_tmp(i,j,kk+KS+1-(KS))
end do
end do
end do
endif
if (including_aerosols) then
do kk = 1,size(totaer_tmp(:,:,:),3)
do j = 1,size(totaer_tmp(:,:,:),2)
do i = 1,size(totaer_tmp(:,:,:),1)
totaer_tmp(i,j,kk) = Optical%totaerooptdep(i,j,kk,7)
end do
end do
end do
do kk = KS,KE
do j = 1,size(tmp2(:,:,:),2)
do i = 1,size(tmp2(:,:,:),1)
tmp2(i,j,kk) = tmp2(i,j,kk) + &
totaer_tmp (i,j,kk+KS+1-(KS))
end do
end do
end do
endif
do kk = KS,KE
do j = 1,size(tt(:,:,:),2)
do i = 1,size(tt(:,:,:),1)
tt(i,j,kk) = EXP(-1.0E+00*tmp2(i,j,kk))* &
Gas_tf%n2o9c (i,j,kk+KS+1-(KS))
end do
end do
end do
endif
!--------------------------------------------------------------------
! calculate or retrieve the source function for the current band.
!--------------------------------------------------------------------
if (n <= 8 + ioffset) then
call looktab (tabsr, ixoe1, dte1, sorc_tmp, KS, KE+1, n)
else
do kk = 1,size(sorc_tmp(:,:,:),3)
do j = 1,size(sorc_tmp(:,:,:),2)
do i = 1,size(sorc_tmp(:,:,:),1)
sorc_tmp(i,j,kk) = sorc(i,j,kk,n-8-ioffset)
end do
end do
end do
endif
!---------------------------------------------------------------------
! retrieve the cfc effect if cfcs are activated.
!---------------------------------------------------------------------
if (Lw_control%do_cfc .and. n >= 9+ioffset) then
call cfc_exact(n-8-ioffset, Optical, cfc_tf)
do kk = KS,KE
do j = 1,size(tt(:,:,:),2)
do i = 1,size(tt(:,:,:),1)
tt(i,j,kk) = tt(i,j,kk)*cfc_tf(i,j,kk)
end do
end do
end do
endif
!------------------------------------------------------------------
! define some near-surface pressure functions that are needed
!------------------------------------------------------------------
do j = 1,size(pfac1(:,:),2)
do i = 1,size(pfac1(:,:),1)
pfac1(i,j) = 0.5*pdfinv2(i,j,KE)*(pflux(i,j,KE+1) - &
press(i,j,KE))*tt(i,j,KE-1)
end do
end do
do j = 1,size(pfac2(:,:),2)
do i = 1,size(pfac2(:,:),1)
pfac2(i,j) = 0.5*pdfinv2(i,j,KE)*(pflux(i,j,KE+1) + &
press(i,j,KE) - 2.0*pflux(i,j,KE))*tt(i,j,KE)
end do
end do
!--------------------------------------------------------------------
! calculate the ground fluxes (?)
!--------------------------------------------------------------------
if (Rad_control%do_totcld_forcing) then
do j = 1,size(fctsgcf(:,:,:),2)
do i = 1,size(fctsgcf(:,:,:),1)
fctsgcf(i,j,n) = tt(i,j,KE)*sorc_tmp(i,j,KE) + &
(pfac1(i,j) + pfac2(i,j))* &
(sorc_tmp(i,j,KE+1) - sorc_tmp(i,j,KE))
end do
end do
do j = 1,size(Lw_diagnostics%fctsg(:,:,:),2)
do i = 1,size(Lw_diagnostics%fctsg(:,:,:),1)
Lw_diagnostics%fctsg(i,j,n) = cldtf(i,j,KE+1,cld_indx_table(n+32-ioffset))* &
fctsgcf(i,j,n)
end do
end do
do j = 1,size(gxctscf(:,:),2)
do i = 1,size(gxctscf(:,:),1)
gxctscf(i,j) = gxctscf(i,j) + fctsgcf(i,j,n)
end do
end do
do j = 1,size(Lw_diagnostics%gxcts(:,:),2)
do i = 1,size(Lw_diagnostics%gxcts(:,:),1)
Lw_diagnostics%gxcts(i,j) = Lw_diagnostics%gxcts(i,j) + &
Lw_diagnostics%fctsg(i,j,n)
end do
end do
else
do j = 1,size(Lw_diagnostics%fctsg(:,:,:),2)
do i = 1,size(Lw_diagnostics%fctsg(:,:,:),1)
Lw_diagnostics%fctsg(i,j,n) = cldtf(i,j,KE+1,cld_indx_table(n+32-ioffset))* &
(tt(i,j,KE)*sorc_tmp(i,j,KE) + &
(pfac1(i,j) + pfac2(i,j))* &
(sorc_tmp(i,j,KE+1) - sorc_tmp(i,j,KE)))
end do
end do
do j = 1,size(Lw_diagnostics%gxcts(:,:),2)
do i = 1,size(Lw_diagnostics%gxcts(:,:),1)
Lw_diagnostics%gxcts(i,j) = Lw_diagnostics%gxcts(i,j) + &
Lw_diagnostics%fctsg(i,j,n)
end do
end do
endif
do j = 1,size(Lw_diagnostics%fctsg(:,:,:),2)
do i = 1,size(Lw_diagnostics%fctsg(:,:,:),1)
Lw_diagnostics%fctsg(i,j,n) = 1.0e-03* &
Lw_diagnostics%fctsg(i,j,n)
end do
end do
!--------------------------------------------------------------------
! include the effect of the cloud transmission function.
!--------------------------------------------------------------------
do kk = KS+1,KE+1
do j = 1,size(ctmp(:,:,:),2)
do i = 1,size(ctmp(:,:,:),1)
ctmp(i,j,kk) = tt(i,j,kk+KS-(KS+1))* &
cldtf(i,j,kk,cld_indx_table(n+32-ioffset))
end do
end do
end do
!---------------------------------------------------------------------
! if diagnostics is on, save each band's contribution separately.
! exctsn is the cool-to-space heating rate for each frequency
! band. fctsg is the "exact" surface flux for each frequency
! band in the 160-1200 cm-1 range.
!---------------------------------------------------------------------
!! the following array only needed when diagnostics on
do kk = KS,KE
do j = 1,size(Lw_diagnostics%exctsn(:,:,:,:),2)
do i = 1,size(Lw_diagnostics%exctsn(:,:,:,:),1)
Lw_diagnostics%exctsn(i,j,kk,n) = sorc_tmp(i,j,kk)* &
(ctmp(i,j,kk+KS+1-(KS)) - ctmp(i,j,kk))
end do
end do
end do
if (Rad_control%do_totcld_forcing) then
do j = 1,size(exctsncf(:,:,:,:),2)
do i = 1,size(exctsncf(:,:,:,:),1)
exctsncf(i,j,KS,n) = sorc_tmp(i,j,KS)* &
(tt(i,j,KS) - 1.0E+00)
end do
end do
do kk = KS+1,KE
do j = 1,size(exctsncf(:,:,:,:),2)
do i = 1,size(exctsncf(:,:,:,:),1)
exctsncf(i,j,kk,n) = sorc_tmp(i,j,kk)* &
(tt(i,j,kk) - tt(i,j,kk+KS-(KS+1)))
end do
end do
end do
endif
!-----------------------------------------------------------------------
! excts is the cool-to-space cooling rate accumulated over
! frequency bands.
!-----------------------------------------------------------------------
do kk = KS,KE
do j = 1,size(Lw_diagnostics%excts(:,:,:),2)
do i = 1,size(Lw_diagnostics%excts(:,:,:),1)
Lw_diagnostics%excts(i,j,kk) = &
Lw_diagnostics%excts(i,j,kk) + &
Lw_diagnostics%exctsn(i,j,kk,n)
end do
end do
end do
if (Rad_control%do_totcld_forcing) then
do j = 1,size(exctscf(:,:,:),2)
do i = 1,size(exctscf(:,:,:),1)
exctscf(i,j,KS) = exctscf(i,j,KS) + &
exctsncf(i,j,KS,n)
end do
end do
do kk = KS+1,KE
do j = 1,size(exctscf(:,:,:),2)
do i = 1,size(exctscf(:,:,:),1)
exctscf(i,j,kk) = exctscf(i,j,kk) + &
exctsncf(i,j,kk,n)
end do
end do
end do
endif
end do
!-----------------------------------------------------------------------
! gxcts is the "exact" surface flux accumulated over
! the frequency bands in the 160-1200 cm-1 range.
! obtain cool-to-space flux at the top by integration of heating
! rates and using cool-to-space flux at the bottom (current value
! of gxcts). note that the pressure quantities and conversion
! factors have not been included either in excts or in gxcts.
! these cancel out, thus reducing computations.
!-----------------------------------------------------------------------
do k=KS,KE
do j = 1,size(Lw_diagnostics%gxcts(:,:),2)
do i = 1,size(Lw_diagnostics%gxcts(:,:),1)
Lw_diagnostics%gxcts(i,j) = Lw_diagnostics%gxcts(i,j) - Lw_diagnostics%excts(i,j,k)
end do
end do
enddo
if (Rad_control%do_totcld_forcing) then
do k=KS,KE
do j = 1,size(gxctscf(:,:),2)
do i = 1,size(gxctscf(:,:),1)
gxctscf(i,j) = gxctscf(i,j) - exctscf(i,j,k)
end do
end do
enddo
endif
!-----------------------------------------------------------------------
! now scale the cooling rate excts by including the pressure
! factor pdfinv and the conversion factor radcon.
!-----------------------------------------------------------------------
!! the following array only needed when diagnostics on
do n=1,NBLY-1
do kk = KS,KE
do j = 1,size(Lw_diagnostics%exctsn(:,:,:,:),2)
do i = 1,size(Lw_diagnostics%exctsn(:,:,:,:),1)
Lw_diagnostics%exctsn(i,j,kk,n) = 1.0e-03*(Lw_diagnostics%exctsn(i,j,kk,n)*radcon_mks* &
pdfinv(i,j,kk))
end do
end do
end do
enddo
if (Rad_control%do_totcld_forcing) then
do n=1,NBLY-1
do kk = KS,KE
do j = 1,size(exctsncf(:,:,:,:),2)
do i = 1,size(exctsncf(:,:,:,:),1)
exctsncf(i,j,kk,n) = 1.0e-03*(exctsncf(i,j,kk,n)*radcon_mks* &
pdfinv(i,j,kk) )
end do
end do
end do
enddo
endif
do kk = KS,KE
do j = 1,size(Lw_diagnostics%excts(:,:,:),2)
do i = 1,size(Lw_diagnostics%excts(:,:,:),1)
Lw_diagnostics%excts(i,j,kk) = (Lw_diagnostics%excts(i,j,kk)*radcon_mks*pdfinv(i,j,kk))*1.0e-03
end do
end do
end do
if (Rad_control%do_totcld_forcing) then
do kk = KS,KE
do j = 1,size(exctscf(:,:,:),2)
do i = 1,size(exctscf(:,:,:),1)
exctscf(i,j,kk) = (exctscf(i,j,kk)*radcon_mks*pdfinv(i,j,kk))*1.0e-03
end do
end do
end do
endif
!--------------------------------------------------------------------
! save the heating rates to be later sent to longwave_driver_mod.
!--------------------------------------------------------------------
do kk = 1,size(cts_sum(:,:,:),3)
do j = 1,size(cts_sum(:,:,:),2)
do i = 1,size(cts_sum(:,:,:),1)
cts_sum(i,j,kk) = Lw_diagnostics%excts(i,j,kk)
end do
end do
end do
do kk = 1,size(cts_sumcf(:,:,:),3)
do j = 1,size(cts_sumcf(:,:,:),2)
do i = 1,size(cts_sumcf(:,:,:),1)
cts_sumcf(i,j,kk) = exctscf(i,j,kk)
end do
end do
end do
!--------------------------------------------------------------------
end subroutine cool_to_space_exact
!####################################################################
!
!
! Subroutine to compute thermal exchange terms and emissivities used
! to obtain the cool-to-space heating rates for all pressure layers.
!
!
! ! E1e290 computes two different quantities.
!
! 1) emissivities used to compute the exchange terms for flux at the
! top of the atmosphere (level KS). (the top layer, isothermal by
! assumption, does not contribute to photon exchanges with other
! layers). these terms are obtained using precomputed e2 functions
! (see ref. (2)).
!
! 2) emissivities used to obtain the cool-to-space heating rates
! for all pressure layers. these are obtained using precomputed
! e1 functions (see ref. (2)).
!
! the frequency ranges for the e2 calculations are 0-560 and 1200-
! 2200 cm-1. the CTS calculations also require calculations in the
! 160-560 cm-1 range. (see refs. (1) and (2)).
!ifdef ch4n2o
!
! if ch4 and n2o are included, the frequency range for emissivities
! is 1400-2200 cm-1, with separate emissivity quantities for the
! 1200-1400 cm-1 range.
!endif ch4n2o
!
! the reason for combining these calculations is that both use
! the same scaled h2o amount (avephi) as input, thus reducing
! some calculation time for obtaining index quantities.
!
! references:
!
! (1) schwarzkopf, m. d., and s. b. fels, "the simplified
! exchange method revisited: an accurate, rapid method for
! computation of infrared cooling rates and fluxes," journal
! of geophysical research, 96 (1981), 9075-9096.
!
! (2) fels, s. b., and m. d. schwarzkopf, "the simplified exchange
! approximation: a new method for radiative transfer
! calculations," journal atmospheric science, 32 (1975),
! 1475-1488.
!
! author: m. d. schwarzkopf
!
! revised: 1/1/93
!
! certified: radiation version 1.0
!
!
! call e1e290 (Atmos_input, e1ctw1, e1ctw2, &
! trans_band1, trans_band2, Optical, tch4n2oe, &
! t4, Lw_diagnostics, cldtf, cld_indx, flx1e1cf, &
! tcfc8)
!
!
! Atmospheric input data to the thermal exchange method
!
!
! Cool to space thermal exchange terms in 0-560 band
!
!
! Cool to space thermal exchange terms in 1200-2200 band
!
!
! transmission functions in band 1
!
!
! transmission function in band 2
!
!
! thermal layer optical path
!
!
! CH4 and N2O transmission functions
!
!
! source function of the top most layer
!
!
! longwave diagnostics variable
!
!
! Cloud transmission function
!
!
! Cloud type index
!
!
! TOA flux due to cloud forcing
!
!
! CFC transmission function (chloroflurocarbons)
!
!
!
subroutine e1e290 (Atmos_input, e1ctw1, e1ctw2, trans_band1, &
trans_band2, Optical, tch4n2oe, t4, Lw_diagnostics, &
cldtf, cld_indx, flx1e1cf, tcfc8, including_aerosols)
!-----------------------------------------------------------------------
!
! E1e290 computes two different quantities.
!
! 1) emissivities used to compute the exchange terms for flux at the
! top of the atmosphere (level KS). (the top layer, isothermal by
! assumption, does not contribute to photon exchanges with other
! layers). these terms are obtained using precomputed e2 functions
! (see ref. (2)).
!
! 2) emissivities used to obtain the cool-to-space heating rates
! for all pressure layers. these are obtained using precomputed
! e1 functions (see ref. (2)).
!
! the frequency ranges for the e2 calculations are 0-560 and 1200-
! 2200 cm-1. the CTS calculations also require calculations in the
! 160-560 cm-1 range. (see refs. (1) and (2)).
!ifdef ch4n2o
!
! if ch4 and n2o are included, the frequency range for emissivities
! is 1400-2200 cm-1, with separate emissivity quantities for the
! 1200-1400 cm-1 range.
!endif ch4n2o
!
! the reason for combining these calculations is that both use
! the same scaled h2o amount (avephi) as input, thus reducing
! some calculation time for obtaining index quantities.
!
! references:
!
! (1) schwarzkopf, m. d., and s. b. fels, "the simplified
! exchange method revisited: an accurate, rapid method for
! computation of infrared cooling rates and fluxes," journal
! of geophysical research, 96 (1981), 9075-9096.
!
! (2) fels, s. b., and m. d. schwarzkopf, "the simplified exchange
! approximation: a new method for radiative transfer
! calculations," journal atmospheric science, 32 (1975),
! 1475-1488.
!
! author: m. d. schwarzkopf
!
! revised: 1/1/93
!
! certified: radiation version 1.0
!-----------------------------------------------------------------------
type(atmos_input_type), intent(in) :: Atmos_input
real, dimension (:,:,:), intent(out) :: e1ctw1, e1ctw2
real, dimension (:,:,:,:), intent(out) :: trans_band1, trans_band2
type(optical_path_type), intent(inout) :: Optical
real, dimension (:,:,:,:), intent(in) :: tch4n2oe
real, dimension(:,:,:), intent(in) :: t4
type(lw_diagnostics_type), intent(inout) :: Lw_diagnostics
real, dimension(:,:,:,:), intent(in) :: cldtf
integer, dimension(:), intent(in) :: cld_indx
real, dimension(:,:), intent(out) :: flx1e1cf
real, dimension (:,:,:), intent(inout) :: tcfc8
logical, intent(in) :: including_aerosols
!---------------------------------------------------------------------
!---------------------------------------------------------------------
! intent(in) variables:
!
! Atmos_input
! tch4n2oe
! t4
! cldtf
! cld_indx
!
! intent(inout) variables:
!
! Optical
! Lw_diagnostics
! tcfc8
!
! intent(out) variables:
!
! e1ctw1
! e1ctw2
! trans_band1
! trans_band2
! flx1e1cf
!
!--------------------------------------------------------------------
!---------------------------------------------------------------------
! local variables
real, dimension (size(trans_band2,1), &
size(trans_band2,2), &
size(trans_band2,3)) :: dte1, dte2,ttmp,ttmp0
integer, dimension (size(trans_band2,1), &
size(trans_band2,2), &
size(trans_band2,3)) :: ixoe1, ixoe2
real, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2), &
size(Atmos_input%temp,3)) :: &
temp, tflux, totaer_tmp, taero8, &
tmp1, tmp2, e1cts1, e1cts2, &
avphilog, dt1, du, dup, du1, dup1
integer, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2), &
size(Atmos_input%temp,3)) :: &
ixo1, iyo, iyop, iyo1, iyop1
real, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2)) :: &
s1a, flxge1, flxge1cf
real, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2), NBTRGE) :: &
flx1e1fcf, flxge1f, flxge1fcf
real, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2), &
size(Atmos_input%temp,3), NBTRGE) :: &
e1cts1f, e1cts2f
integer :: i,j,k,kk,m
!---------------------------------------------------------------------
! local variables
!
! dte1
! dte2
! ixoe1
! ixoe2
! temp
! tflux
! totaer_tmp
! taero8
! tmp1
! tmp2
! e1cts1
! e1cts2
! avphilog
! dt1
! du
! dup
! ixo1
! iyo
! iyop
! s1a
! flxge1
! flxge1cf
! flx1e1fcf
! flxge1f
! flxge1fcf
! e1cts1f
! e1cts2f
! k,m
!
!----------------------------------------------------------------------
do kk = 1,size(tflux(:,:,:),3)
do j = 1,size(tflux(:,:,:),2)
do i = 1,size(tflux(:,:,:),1)
tflux(i,j,kk) = Atmos_input%tflux(i,j,kk)
end do
end do
end do
do kk = 1,size(temp(:,:,:),3)
do j = 1,size(temp(:,:,:),2)
do i = 1,size(temp(:,:,:),1)
temp(i,j,kk) = Atmos_input%temp(i,j,kk)
end do
end do
end do
!---------------------------------------------------------------------
! obtain the "exchange" emissivities as a function of temperature
! (fxo) and water amount (fyo). the temperature indices have
! been obtained in longwave_setup_mod.
!-------------------------------------------------------------------
call locate_in_table(temp_1, temp, dte1, ixoe1, KS, KE+1)
call locate_in_table(temp_1, tflux, dte2, ixoe2, KS, KE+1)
do kk = KS,KE
do j = 1,size(ixoe2(:,:,:),2)
do i = 1,size(ixoe2(:,:,:),1)
ixoe2(i,j,kk) = ixoe2(i,j,kk+KS+1-(KS))
end do
end do
end do
do kk = KS,KE
do j = 1,size(dte2(:,:,:),2)
do i = 1,size(dte2(:,:,:),1)
dte2(i,j,kk) = dte2 (i,j,kk+KS+1-(KS))
end do
end do
end do
do j = 1,size(ixoe2(:,:,:),2)
do i = 1,size(ixoe2(:,:,:),1)
ixoe2(i,j,KE+1) = ixoe1(i,j,KE)
end do
end do
do j = 1,size(dte2(:,:,:),2)
do i = 1,size(dte2(:,:,:),1)
dte2(i,j,KE+1) = dte1 (i,j,KE)
end do
end do
if (Lw_control%do_h2o) then
do kk = KS,KE+1
do j = 1,size(avphilog(:,:,:),2)
do i = 1,size(avphilog(:,:,:),1)
avphilog(i,j,kk) = LOG10(Optical%avephi(i,j,kk))
end do
end do
end do
call locate_in_table (mass_1, avphilog, du, iyo, KS, KE+1)
iyo(:,:,KS:KE+1) = iyo(:,:,KS:KE+1) + 1
call looktab (tab2, ixoe2, iyo, dte2, du, &
trans_band2(:,:,:,1), KS, KE+1)
else
iyo1(:,:,KS:KE+1) = 1
du1(:,:,KS:KE+1) = 0.0
call looktab (tab2, ixoe2, iyo1, dte2, du1, &
trans_band2(:,:,:,1), KS, KE+1)
endif
!-----------------------------------------------------------------------
! the special case emiss(:,:,KE+1) for layer KE+1 is obtained by
! averaging the values for KE and KE+1.
!---------------------------------------------------------------------
do j = 1,size(trans_band2(:,:,:,:),2)
do i = 1,size(trans_band2(:,:,:,:),1)
trans_band2(i,j,KE+1,1) = 0.5E+00*(trans_band2(i,j,KE,1) + &
trans_band2(i,j,KE+1, 1))
end do
end do
ttmp(:,:,:) = trans_band2(:,:,:,1)
do kk = KS+1,KE
do j = 1,size(trans_band2(:,:,:,:),2)
do i = 1,size(trans_band2(:,:,:,:),1)
trans_band2(i,j,kk,1) = ttmp(i,j,kk-1)
end do
end do
end do
!---------------------------------------------------------------------
! perform calculations for the e1 function. the terms involving top
! layer du are not known. we use index two to represent index one
! in previous calculations. (now 3)
!--------------------------------------------------------------------
do j = 1,size(iyop(:,:,:),2)
do i = 1,size(iyop(:,:,:),1)
iyop(i,j,KS) = 2
end do
end do
do kk = KS+1,KE+1
do j = 1,size(iyop(:,:,:),2)
do i = 1,size(iyop(:,:,:),1)
iyop(i,j,kk) = iyo(i,j,kk+KS-(KS+1))
end do
end do
end do
do j = 1,size(dup(:,:,:),2)
do i = 1,size(dup(:,:,:),1)
dup(i,j,KS) = 0.0E+00
end do
end do
do kk = KS+1,KE+1
do j = 1,size(dup(:,:,:),2)
do i = 1,size(dup(:,:,:),1)
dup(i,j,kk) = du (i,j,kk+KS-(KS+1))
end do
end do
end do
do k=KS,KE+1
do j = 1,size(ixo1(:,:,:),2)
do i = 1,size(ixo1(:,:,:),1)
ixo1(i,j,k) = ixoe1(i,j,KS)
end do
end do
do j = 1,size(dt1(:,:,:),2)
do i = 1,size(dt1(:,:,:),1)
dt1(i,j,k) = dte1 (i,j,KS)
end do
end do
enddo
if (Lw_control%do_h2o) then
!-----------------------------------------------------------------------
! e1flx(:,:,KS) equals e1cts1(:,:,KS).
!-----------------------------------------------------------------------
call looktab (tab1, ixoe1, iyop, dte1, dup, e1cts1, KS, KE+1)
call looktab (tab1, ixoe1, iyo, dte1, du, e1cts2, KS, KE)
call looktab (tab1, ixo1, iyop, dt1, dup, &
trans_band1(:,:,:,1), KS, KE+1)
call looktab (tab1w, ixoe1, iyop, dte1, dup, e1ctw1, KS, KE+1)
call looktab (tab1w, ixoe1, iyo, dte1, du, e1ctw2, KS, KE)
else
iyop1(:,:,:) = 1
dup1 (:,:,:) = 0.0
call looktab (tab1, ixoe1, iyop1, dte1, dup1, e1cts1, KS, KE+1)
call looktab (tab1, ixoe1, iyo1, dte1, du1, e1cts2, KS, KE)
call looktab (tab1, ixo1, iyop1, dt1, dup1, &
trans_band1(:,:,:,1), KS, KE+1)
call looktab (tab1w, ixoe1, iyop1, dte1, dup1, e1ctw1, KS, KE+1)
call looktab (tab1w, ixoe1, iyo1, dte1, du1, e1ctw2, KS, KE)
endif
!--------------------------------------------------------------------
! calculations with ch4 and n2o require NBTRGE separate emissivity
! bands for h2o.
!--------------------------------------------------------------------
if (NBTRGE > 0) then
do m=1,NBTRGE
if (Lw_control%do_h2o) then
do kk = KS,KE+1
do j = 1,size(avphilog(:,:,:),2)
do i = 1,size(avphilog(:,:,:),1)
avphilog(i,j,kk) = LOG10(Optical%avephif(i,j,kk,m))
end do
end do
end do
call locate_in_table (mass_1, avphilog, du, iyo, KS , KE+1)
iyo(:,:,KS:KE+1) = iyo(:,:,KS:KE+1) + 1
do j = 1,size(iyop(:,:,:),2)
do i = 1,size(iyop(:,:,:),1)
iyop(i,j,KS) = 2
end do
end do
do kk = KS+1,KE+1
do j = 1,size(iyop(:,:,:),2)
do i = 1,size(iyop(:,:,:),1)
iyop(i,j,kk) = iyo(i,j,kk+KS-(KS+1))
end do
end do
end do
do j = 1,size(dup(:,:,:),2)
do i = 1,size(dup(:,:,:),1)
dup(i,j,KS) = 0.0E+00
end do
end do
do kk = KS+1,KE+1
do j = 1,size(dup(:,:,:),2)
do i = 1,size(dup(:,:,:),1)
dup(i,j,kk) = du (i,j,kk+KS-(KS+1))
end do
end do
end do
call looktab (tab2a, ixoe2, iyo, dte2, du, &
trans_band2(:,:,:,6+m), KS, KE+1, m)
call looktab (tab1a, ixoe1, iyop, dte1, dup, &
e1cts1f(:,:,:,m), KS, KE+1, m)
call looktab (tab1a, ixoe1, iyo, dte1, du, e1cts2f(:,:,:,m),&
KS, KE, m)
call looktab (tab1a, ixo1, iyop, dt1, dup, &
trans_band1(:,:,:,6+m), KS, KE+1, m)
else
iyo1(:,:,KS:KE+1) = 1
iyop1(:,:,KS:KE+1) = 1
du1(:,:,KS:KE+1) = 0.0
dup1(:,:,KS:KE+1) = 0.0
call looktab (tab2a, ixoe2, iyo1, dte2, du1, &
trans_band2(:,:,:,6+m), KS, KE+1, m)
call looktab (tab1a, ixoe1, iyop1, dte1, dup1, &
e1cts1f(:,:,:,m), KS, KE+1, m)
call looktab (tab1a, ixoe1, iyo1, dte1, du1, &
e1cts2f(:,:,:,m), KS, KE, m)
call looktab (tab1a, ixo1, iyop1, dt1, dup1, &
trans_band1(:,:,:,6+m), KS, KE+1, m)
endif
enddo
!--------------------------------------------------------------------
! the special case emissf(:,:,KE+1,m) for layer KE+1 is obtained by
! averaging the values for KE and KE+1.
!--------------------------------------------------------------------
do m=1,NBTRGE
ttmp(:,:,:) = trans_band2(:,:,:,6+m)
ttmp0(:,:,:) = trans_band2(:,:,:,6+m)
do j = 1,size(trans_band2(:,:,:,:),2)
do i = 1,size(trans_band2(:,:,:,:),1)
trans_band2(i,j,KE+1,6+m) = 0.5E+00* &
(ttmp(i,j,KE) + &
ttmp0(i,j,KE+1))
end do
end do
ttmp(:,:,:)=trans_band2(:,:,:,6+m)
do kk = KS+1,KE
do j = 1,size(trans_band2(:,:,:,:),2)
do i = 1,size(trans_band2(:,:,:,:),1)
trans_band2(i,j,kk,6+m) = ttmp(i,j,kk+KS-(KS+1))
end do
end do
end do
enddo
endif
!---------------------------------------------------------------------
! add the effects of other radiative gases on these flux arrays.
! the lbl transmissivities imply (at present) NBTRG = NBTRGE = 1).
! thus, tch4e and tn2oe are obtained directly from the transmission
! functions.
!----------------------------------------------------------------------
if (NBTRGE > 0) then
if (Lw_control%do_ch4 .or. Lw_control%do_n2o) then
do kk = KS+1,KE+1
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kk,6+1) = trans_band1(i,j,kk,6+1)* &
tch4n2oe(i,j,kk,1)
end do
end do
end do
do kk = KS,KE+1
do j = 1,size(e1cts1f(:,:,:,:),2)
do i = 1,size(e1cts1f(:,:,:,:),1)
e1cts1f(i,j,kk,1) = e1cts1f(i,j,kk,1)* &
tch4n2oe(i,j,kk,1)
end do
end do
end do
do kk = KS,KE
do j = 1,size(e1cts2f(:,:,:,:),2)
do i = 1,size(e1cts2f(:,:,:,:),1)
e1cts2f(i,j,kk, 1) = e1cts2f(i,j,kk,1)* &
tch4n2oe(i,j,kk+KS+1-(KS),1)
end do
end do
end do
do kk = KS+1,KE+1
do j = 1,size(trans_band2(:,:,:,:),2)
do i = 1,size(trans_band2(:,:,:,:),1)
trans_band2(i,j,kk, 6+1) = trans_band2(i,j,kk,6+1)* &
tch4n2oe(i,j,kk,1)
end do
end do
end do
endif
!----------------------------------------------------------------------
! add cfc transmissivities if species which absorb in this fre-
! quency range ( presently index 8) are present.
!----------------------------------------------------------------------
!----------------------------------------------------------------------
if (Lw_control%do_cfc) then
call cfc_indx8 (8, Optical, tcfc8)
do kk = KS+1,KE+1
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kk,6+1) = trans_band1(i,j,kk,6+1)* &
tcfc8(i,j,kk)
end do
end do
end do
do kk = KS,KE+1
do j = 1,size(e1cts1f(:,:,:,:),2)
do i = 1,size(e1cts1f(:,:,:,:),1)
e1cts1f(i,j,kk,1) = e1cts1f(i,j,kk,1)* &
tcfc8(i,j,kk)
end do
end do
end do
do kk = KS,KE
do j = 1,size(e1cts2f(:,:,:,:),2)
do i = 1,size(e1cts2f(:,:,:,:),1)
e1cts2f(i,j,kk, 1) = e1cts2f(i,j,kk,1)* &
tcfc8(i,j,kk+KS+1-(KS))
end do
end do
end do
do kk = KS+1,KE+1
do j = 1,size(trans_band2(:,:,:,:),2)
do i = 1,size(trans_band2(:,:,:,:),1)
trans_band2(i,j,kk,6+1) = trans_band2(i,j,kk,6+1)* &
tcfc8(i,j,kk)
end do
end do
end do
endif
!----------------------------------------------------------------------
! compute aerosol transmission function for 1200-1400 cm-1 region
!----------------------------------------------------------------------
if (including_aerosols) then
do kk = 1,size(totaer_tmp(:,:,:),3)
do j = 1,size(totaer_tmp(:,:,:),2)
do i = 1,size(totaer_tmp(:,:,:),1)
totaer_tmp(i,j,kk) = Optical%totaerooptdep(i,j,kk,8)
end do
end do
end do
do kk = KS,KE+1
do j = 1,size(taero8(:,:,:),2)
do i = 1,size(taero8(:,:,:),1)
taero8(i,j,kk) = EXP(-1.0E+00*totaer_tmp(i,j,kk))
end do
end do
end do
do kk = KS+1,KE+1
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kk,6+1) = trans_band1(i,j,kk,6+1)* &
taero8(i,j,kk)
end do
end do
end do
do kk = KS,KE+1
do j = 1,size(e1cts1f(:,:,:,:),2)
do i = 1,size(e1cts1f(:,:,:,:),1)
e1cts1f(i,j,kk,1) = e1cts1f(i,j,kk,1)* &
taero8(i,j,kk)
end do
end do
end do
do kk = KS,KE
do j = 1,size(e1cts2f(:,:,:,:),2)
do i = 1,size(e1cts2f(:,:,:,:),1)
e1cts2f(i,j,kk, 1) = e1cts2f(i,j,kk,1)* &
taero8(i,j,kk+KS+1-(KS))
end do
end do
end do
do kk = KS+1,KE+1
do j = 1,size(trans_band2(:,:,:,:),2)
do i = 1,size(trans_band2(:,:,:,:),1)
trans_band2(i,j,kk,6+1) = trans_band2(i,j,kk,6+1)* &
taero8(i,j,kk)
end do
end do
end do
endif
endif ! (NBTRGE > 0)
!-----------------------------------------------------------------------
! obtain the flux at the top of the atmosphere in the 0-160,
! 1200-2200 cm-1 frequency ranges, where heating rates and fluxes
! are derived from h2o emissivity calculations (flx1e1) by:
! 1) obtaining the surface flux (flxge1); 2) summing the
! emissivity flux divergence for these ranges (tmp1) over all
! pressure layers.
!#ifdef ch4n2o
! if the 1200-1400 cm-1 range is computed separately, flux calcu-
! lations are done separately in this range, then combined with
! those from the other frequency range.
!#endif ch4n2o
!----------------------------------------------------------------------
do j = 1,size(s1a(:,:),2)
do i = 1,size(s1a(:,:),1)
s1a(i,j) = t4(i,j,KE+1)*(e1cts1(i,j,KE+1) - e1ctw1(i,j,KE+1))
end do
end do
do j = 1,size(flxge1(:,:),2)
do i = 1,size(flxge1(:,:),1)
flxge1(i,j) = s1a(i,j)*cldtf(i,j,KE+1,1)
end do
end do
do kk = KS,KE
do j = 1,size(tmp1(:,:,:),2)
do i = 1,size(tmp1(:,:,:),1)
tmp1(i,j,kk) = t4(i,j,kk)* &
(e1cts1(i,j,kk) - e1ctw1(i,j,kk))
end do
end do
end do
do kk = KS,KE
do j = 1,size(tmp2(:,:,:),2)
do i = 1,size(tmp2(:,:,:),1)
tmp2(i,j,kk) = t4(i,j,kk)* &
(e1cts2(i,j,kk) - e1ctw2(i,j,kk))
end do
end do
end do
do j = 1,size(Lw_diagnostics%flx1e1(:,:),2)
do i = 1,size(Lw_diagnostics%flx1e1(:,:),1)
Lw_diagnostics%flx1e1(i,j) = flxge1(i,j)
end do
end do
do k=KS,KE
do j = 1,size(Lw_diagnostics%flx1e1(:,:),2)
do i = 1,size(Lw_diagnostics%flx1e1(:,:),1)
Lw_diagnostics%flx1e1(i,j) = Lw_diagnostics%flx1e1(i,j) + tmp1(i,j,k)*cldtf(i,j,k,1) - &
tmp2(i,j,k)*cldtf(i,j,k+1,1)
end do
end do
enddo
if (Rad_control%do_totcld_forcing) then
do j = 1,size(flxge1cf(:,:),2)
do i = 1,size(flxge1cf(:,:),1)
flxge1cf(i,j) = s1a(i,j)
end do
end do
do j = 1,size(flx1e1cf(:,:),2)
do i = 1,size(flx1e1cf(:,:),1)
flx1e1cf(i,j) = flxge1cf(i,j)
end do
end do
do k=KS,KE
do j = 1,size(flx1e1cf(:,:),2)
do i = 1,size(flx1e1cf(:,:),1)
flx1e1cf(i,j) = flx1e1cf(i,j) + tmp1(i,j,k) - tmp2(i,j,k)
end do
end do
enddo
endif
if (NBTRGE > 0) then
do m=1,NBTRGE
do j = 1,size(s1a(:,:),2)
do i = 1,size(s1a(:,:),1)
s1a(i,j) = t4(i,j,KE+1)*e1cts1f(i,j,KE+1,m)
end do
end do
do j = 1,size(flxge1f(:,:,:),2)
do i = 1,size(flxge1f(:,:,:),1)
flxge1f(i,j,m) = s1a(i,j)*cldtf(i,j,KE+1,cld_indx(7))
end do
end do
do j = 1,size(Lw_diagnostics%flx1e1f(:,:,:),2)
do i = 1,size(Lw_diagnostics%flx1e1f(:,:,:),1)
Lw_diagnostics%flx1e1f(i,j,m) = flxge1f(i,j,m)
end do
end do
do k=KS,KE
do j = 1,size(tmp1(:,:,:),2)
do i = 1,size(tmp1(:,:,:),1)
tmp1(i,j,k) = t4(i,j,k)*e1cts1f(i,j,k,m)
end do
end do
do j = 1,size(tmp2(:,:,:),2)
do i = 1,size(tmp2(:,:,:),1)
tmp2(i,j,k) = t4(i,j,k)*e1cts2f(i,j,k,m)
end do
end do
do j = 1,size(Lw_diagnostics%flx1e1f(:,:,:),2)
do i = 1,size(Lw_diagnostics%flx1e1f(:,:,:),1)
Lw_diagnostics%flx1e1f(i,j,m) = &
Lw_diagnostics%flx1e1f(i,j,m) + tmp1(i,j,k)* &
cldtf(i,j,k,cld_indx(7)) - tmp2(i,j,k)* &
cldtf(i,j,k+1,cld_indx(7))
end do
end do
end do
end do
do m=1,NBTRGE
do j = 1,size(Lw_diagnostics%flx1e1(:,:),2)
do i = 1,size(Lw_diagnostics%flx1e1(:,:),1)
Lw_diagnostics%flx1e1(i,j) = &
Lw_diagnostics%flx1e1(i,j) + &
Lw_diagnostics%flx1e1f(i,j,m)
end do
end do
enddo
if (Rad_control%do_totcld_forcing) then
do m=1,NBTRGE
do j = 1,size(flxge1fcf(:,:,:),2)
do i = 1,size(flxge1fcf(:,:,:),1)
flxge1fcf(i,j,m) = s1a(i,j)
end do
end do
do j = 1,size(flx1e1fcf(:,:,:),2)
do i = 1,size(flx1e1fcf(:,:,:),1)
flx1e1fcf(i,j,m) = s1a(i,j)
end do
end do
do k=KS,KE
do j = 1,size(flx1e1fcf(:,:,:),2)
do i = 1,size(flx1e1fcf(:,:,:),1)
flx1e1fcf(i,j,m) = flx1e1fcf(i,j,m) + tmp1(i,j,k) - &
tmp2(i,j,k)
end do
end do
end do
end do
do m=1,NBTRGE
do j = 1,size(flx1e1cf(:,:),2)
do i = 1,size(flx1e1cf(:,:),1)
flx1e1cf(i,j) = flx1e1cf(i,j) + flx1e1fcf(i,j,m)
end do
end do
enddo
endif
endif ! (ntrge > 0)
!----------------------------------------------------------------------
end subroutine e1e290
!###################################################################
!
!
! e290 computes the exchange terms in the flux equation for longwave
! radiation for all terms except the exchange with the top of the
! atmosphere.
!
!
! e290 computes the exchange terms in the flux equation for longwave
! radiation for all terms except the exchange with the top of the
! atmosphere. the method is a table lookup on a pre-computed e2
! function (defined in reference (2)). calculation are done in the
! frequency range: 0-560, 1200-2200 cm-1 for q(approximate).
! motivation for these calculations is in references (1) and (2).
!
! references:
!
! (1) schwarzkopf, m. d., and s. b. fels, "the simplified
! exchange method revisited: an accurate, rapid method for
! computation of infrared cooling rates and fluxes," journal
! of geophysical research, 96 (1981), 9075-9096.
!
! (2) fels, s. b., and m. d. schwarzkopf, "the simplified exchange
! approximation: a new method for radiative transfer
! calculations," journal atmospheric science, 32 (1975),
! 1475-1488.
!
!
! call e290 (Atmos_input, k, trans_band2, trans_band1, Optical, tch4n2oe, &
! tcfc8)
!
!
! Atmospheric input data
!
!
! Starting vertical level k to compute exchange terms
!
!
! Transmission funciton in band 1200-2200
!
!
! Transmission function in band 0-560
!
!
! Optical depth of the thermal layers
!
!
! CH4 and N2O transmission function
!
!
! CFC transmission function
!
!
!
subroutine e290 (Atmos_input, k, trans_band2, trans_band1, Optical, &
tch4n2oe, tcfc8, including_aerosols)
!-----------------------------------------------------------------------
!
! e290 computes the exchange terms in the flux equation for longwave
! radiation for all terms except the exchange with the top of the
! atmosphere. the method is a table lookup on a pre-computed e2
! function (defined in reference (2)). calculation are done in the
! frequency range: 0-560, 1200-2200 cm-1 for q(approximate).
! motivation for these calculations is in references (1) and (2).
!
! references:
!
! (1) schwarzkopf, m. d., and s. b. fels, "the simplified
! exchange method revisited: an accurate, rapid method for
! computation of infrared cooling rates and fluxes," journal
! of geophysical research, 96 (1981), 9075-9096.
!
! (2) fels, s. b., and m. d. schwarzkopf, "the simplified exchange
! approximation: a new method for radiative transfer
! calculations," journal atmospheric science, 32 (1975),
! 1475-1488.
!
! author: c. h. goldberg
!
! revised: 1/1/93
!
! certified: radiation version 1.0
!
!---------------------------------------------------------------------
type(atmos_input_type), intent(in) :: Atmos_input
integer, intent(in) :: k
real, dimension(:,:,:,:), intent(inout) :: trans_band1, trans_band2
type(optical_path_type), intent(inout) :: Optical
real, dimension(:,:,:,:), intent(in) :: tch4n2oe
real, dimension(:,:,:), intent(inout) :: tcfc8
logical, intent(in) :: including_aerosols
!----------------------------------------------------------------------
!-------------------------------------------------------------------
! intent(in) variables:
!
! Atmos_input
! k
! tch4n2oe
!
! intent(inout) variables:
!
! trans_band1
! trans_band2
! Optical
! tcfc8
!
!---------------------------------------------------------------------
!-------------------------------------------------------------------
! local variables
!-------------------------------------------------------------------
integer :: i,j,kk,kp, m
real, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2), &
size(Atmos_input%temp,3)) :: &
temp, tflux, totaer_tmp, taero8, &
taero8kp, avphilog, dtk, du, du1
integer, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2), &
size(Atmos_input%temp,3)) :: &
ixok, iyo, iyo1
real, dimension (size(trans_band2,1), &
size(trans_band2,2), &
size(trans_band2,3)) :: dte1, dte2
integer, dimension (size(trans_band2,1), &
size(trans_band2,2), &
size(trans_band2,3)) :: ixoe1, ixoe2
real, dimension (size(trans_band2,1), &
size(trans_band2,2), &
size(trans_band2,3)) ::ttmp
real, dimension (size(trans_band1,1), &
size(trans_band1,2), &
size(trans_band1,3)) :: ttmp0
!-------------------------------------------------------------------
! local variables:
!
! temp
! tflux
! totaer_tmp
! taero8
! taero8kp
! avphilog
! dtk
! du
! du1
! ixok
! iyo
! iyo1
! dte1
! dte2
! ixoe1
! ixoe2
! kp,m
!
!---------------------------------------------------------------------
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
temp = Atmos_input%temp
tflux = Atmos_input%tflux
!-----------------------------------------------------------------------
! obtain the "exchange" emissivities as a function of temperature
! (fxo) and water amount (avephi). the temperature indices have
! been obtained in Lwrad. calculations are for flux level k, with
! kp = k+1 to KE+1. the case k=KS is excluded (done in E1e290).
! calculations are also made for flux levels k to KE, for
! contributions from flux level k. in this case, the temperature
! index (ixok) represents tflux(:,:,k-1); the water index (iyo)
! has the same values as in the case with varying kp.
!---------------------------------------------------------------------
!----------------------------------------------------------------------
call locate_in_table(temp_1, temp, dte1, ixoe1, KS, KE+1)
call locate_in_table(temp_1, tflux, dte2, ixoe2, KS, KE+1)
ttmp(:,:,:) = ixoe2(:,:,:)
do kk = KS,KE
do j = 1,size(ixoe2(:,:,:),2)
do i = 1,size(ixoe2(:,:,:),1)
ixoe2(i,j,kk) = ttmp(i,j,kk+KS+1-(KS))
end do
end do
end do
ttmp(:,:,:)=dte2 (:,:,:)
do kk = KS,KE
do j = 1,size(dte2(:,:,:),2)
do i = 1,size(dte2(:,:,:),1)
dte2(i,j,kk) = ttmp(i,j,kk+KS+1-(KS))
end do
end do
end do
do j = 1,size(ixoe2(:,:,:),2)
do i = 1,size(ixoe2(:,:,:),1)
ixoe2(i,j,KE+1) = ixoe1(i,j,KE)
end do
end do
do j = 1,size(dte2(:,:,:),2)
do i = 1,size(dte2(:,:,:),1)
dte2(i,j,KE+1) = dte1 (i,j,KE)
end do
end do
do kp=k,KE
do j = 1,size(ixok(:,:,:),2)
do i = 1,size(ixok(:,:,:),1)
ixok(i,j,kp) = ixoe2(i,j,k-1)
end do
end do
do j = 1,size(dtk(:,:,:),2)
do i = 1,size(dtk(:,:,:),1)
dtk(i,j,kp) = dte2 (i,j,k-1)
end do
end do
end do
if (Lw_control%do_h2o) then
do kk = k,KE+1
do j = 1,size(avphilog(:,:,:),2)
do i = 1,size(avphilog(:,:,:),1)
avphilog(i,j,kk) = LOG10(Optical%avephi(i,j,kk))
end do
end do
end do
call locate_in_table (mass_1, avphilog, du, iyo,k, KE+1)
iyo(:,:,k:KE+1) = iyo(:,:,k:KE+1) + 1
call looktab (tab2, ixoe2, iyo, dte2, du, &
trans_band2(:,:,:,1), k, KE+1)
call looktab (tab2, ixok, iyo, dtk, du, &
trans_band1(:,:,:,1), k, KE)
else
iyo1(:,:,k:KE+1) = 1
du1(:,:,k:KE+1) = 0.0
call looktab (tab2, ixoe2, iyo1, dte2, du1, &
trans_band2(:,:,:,1), k, KE+1)
call looktab (tab2, ixok, iyo1, dtk, du1, &
trans_band1(:,:,:,1), k, KE)
endif
ttmp0(:,:,:)=trans_band1(:,:,:,1)
do kk = k+1,KE+1
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kk,1) = ttmp0(i,j,kk+k-(k+1))
end do
end do
end do
!--------------------------------------------------------------------
! the special case emiss(:,:,KE) for layer KE is obtained by
! averaging the values for KE and KE+1. note that emiss(:,:,KE+1)
! is not useful after this point.
!-------------------------------------------------------------------
do j = 1,size(trans_band2(:,:,:,:),2)
do i = 1,size(trans_band2(:,:,:,:),1)
trans_band2(i,j,KE+1,1) = 0.5E+00*(trans_band2(i,j,KE,1) + &
trans_band2(i,j,KE+1,1))
end do
end do
ttmp(:,:,:)=trans_band2(:,:,:,1)
do kk = k+1,KE
do j = 1,size(trans_band2(:,:,:,:),2)
do i = 1,size(trans_band2(:,:,:,:),1)
trans_band2(i,j,kk,1) = ttmp(i,j,kk+k-(k+1))
end do
end do
end do
!--------------------------------------------------------------------
! calculations with ch4 and n2o require NBTRGE separate emissivity
! bands for h2o. reqults are in emissf (flux level k) and
! emissbf (other levels).
!-------------------------------------------------------------------
if (nbtrge > 0) then
do m=1,NBTRGE
if (Lw_control%do_h2o) then
do kk = k,KE+1
do j = 1,size(avphilog(:,:,:),2)
do i = 1,size(avphilog(:,:,:),1)
avphilog(i,j,kk) = LOG10(Optical%avephif(i,j,kk,m))
end do
end do
end do
call locate_in_table (mass_1, avphilog, du, iyo, k, KE+1)
iyo(:,:,k:KE+1) = iyo(:,:,k:KE+1) + 1
call looktab (tab2a, ixoe2, iyo, dte2, du, &
trans_band2(:,:,:,6+m), k, KE+1, m)
call looktab (tab2a, ixok, iyo, dtk, du, &
trans_band1(:,:,:,6+m), k, KE, m)
else
iyo1(:,:,k:KE+1) = 1
du1(:,:,k:KE+1) = 0.0
call looktab (tab2a, ixoe2, iyo1, dte2, du1, &
trans_band2(:,:,:,6+m), k, KE+1, m)
call looktab (tab2a, ixok, iyo1, dtk, du1, &
trans_band1(:,:,:,6+m), k, KE, m)
endif
ttmp0(:,:,:)=trans_band1(:,:,:,6+m)
do kk = k+1,KE+1
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kk,6+m) = ttmp0(i,j,kk+k-(k+1))
end do
end do
end do
enddo
!----------------------------------------------------------------------
! the special case emissf(:,:,KE) for layer KE is obtained by
! averaging the values for KE and KE+1. note that emissf(:,:,KE+1,m)
! is not useful after this point.
!----------------------------------------------------------------------
do m=1,NBTRGE
do j = 1,size(trans_band2(:,:,:,:),2)
do i = 1,size(trans_band2(:,:,:,:),1)
trans_band2(i,j,KE+1,6+m) = 0.5E+00* &
(trans_band2(i,j,KE,6+m) + &
trans_band2(i,j,KE+1,6+m))
end do
end do
ttmp(:,:,:)=trans_band2(:,:,:,6+m)
do kk = k+1,KE
do j = 1,size(trans_band2(:,:,:,:),2)
do i = 1,size(trans_band2(:,:,:,:),1)
trans_band2(i,j,kk,6+m) = ttmp(i,j,kk+k-(k+1))
end do
end do
end do
enddo
endif
!----------------------------------------------------------------------
! add the effects of other radiative gases on these flux arrays.
! the lbl transmissivities imply (at present) NBTRG = NBTRGE = 1).
! thus, tch4e and tn2oe are obtained directly from the transmission
! functions.
!---------------------------------------------------------------------
if (nbtrge > 0) then
if (Lw_control%do_ch4 .or. Lw_control%do_n2o) then
do kk = k+1,KE+1
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kk,6+1) = trans_band1(i,j,kk,6+1)* &
tch4n2oe(i,j,kk,1)
end do
end do
end do
do kk = k+1,KE+1
do j = 1,size(trans_band2(:,:,:,:),2)
do i = 1,size(trans_band2(:,:,:,:),1)
trans_band2(i,j,kk,6+1) = trans_band2(i,j,kk,6+1)*tch4n2oe(i,j,kk,1)
end do
end do
end do
endif
!--------------------------------------------------------------------
! add cfc transmissivities if species which absorb in this fre-
! quency range are present.
!--------------------------------------------------------------------
if (Lw_control%do_cfc) then
call cfc_indx8_part (8, Optical, tcfc8, k)
do kk = k+1,KE+1
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kk,6+1) = trans_band1(i,j,kk,6+1)*tcfc8(i,j,kk)
end do
end do
end do
do kk = k+1,KE+1
do j = 1,size(trans_band2(:,:,:,:),2)
do i = 1,size(trans_band2(:,:,:,:),1)
trans_band2(i,j,kk,6+1) = trans_band2(i,j,kk,6+1)*tcfc8(i,j,kk)
end do
end do
end do
endif
!--------------------------------------------------------------------
! compute aerosol transmission function for 1200-1400 cm-1 region
! (as quotient of 2 exponentials)
! taero8kp(k) contains the (k+1,k) transmissivities for all k
! in the 1200-1400 cm-1 frequency range.
!---------------------------------------------------------------------
if (including_aerosols) then
do kk = 1,size(totaer_tmp(:,:,:),3)
do j = 1,size(totaer_tmp(:,:,:),2)
do i = 1,size(totaer_tmp(:,:,:),1)
totaer_tmp(i,j,kk) = Optical%totaerooptdep(i,j,kk,8)
end do
end do
end do
do kk = KS,KE+1
do j = 1,size(taero8(:,:,:),2)
do i = 1,size(taero8(:,:,:),1)
taero8(i,j,kk) = EXP(-1.0E+00*totaer_tmp(i,j,kk))
end do
end do
end do
do kp = k+1,KE+1
do j = 1,size(taero8kp(:,:,:),2)
do i = 1,size(taero8kp(:,:,:),1)
taero8kp(i,j,kp) = taero8(i,j,kp)/taero8(i,j,k)
end do
end do
enddo
do kk = k+1,KE+1
do j = 1,size(trans_band1(:,:,:,:),2)
do i = 1,size(trans_band1(:,:,:,:),1)
trans_band1(i,j,kk,6+1) = trans_band1(i,j,kk,6+1)* &
taero8kp(i,j,kk)
end do
end do
end do
do kk = k+1,KE+1
do j = 1,size(trans_band2(:,:,:,:),2)
do i = 1,size(trans_band2(:,:,:,:),1)
trans_band2(i,j,kk,6+1) = trans_band2(i,j,kk,6+1)* &
taero8kp(i,j,kk)
end do
end do
end do
endif
endif ! (nbtrge > 0)
!--------------------------------------------------------------------
end subroutine e290
!####################################################################
!
!
! Subroutine to compute thermal layer emissivity using pre computed
! look up tables
!
!
! call esfc (Atmos_input, emspec, Optical, &
! emspecf, tch4n2oe, tcfc8 )
!
!
! Atmospheric input data such as temperature and flux level temp
!
!
! Emissivity of thermal layers
!
!
! Optical depth of thermal layers
!
!
! Emissivity of thermal layers including effects of minor gas species
!
!
! CH4 and N2O transmission function
!
!
! CFC transmission function
!
!
!
subroutine esfc (Atmos_input, emspec, Optical, emspecf, tch4n2oe, tcfc8 )
!--------------------------------------------------------------------
!
!--------------------------------------------------------------------
type(atmos_input_type), intent(in) :: Atmos_input
real, dimension (:,:,:), intent(out) :: emspec
type(optical_path_type), intent(inout) :: Optical
real, dimension (:,:,:,:), intent(out) :: emspecf
real, dimension (:,:,:,:), intent(in) :: tch4n2oe
real, dimension (:,:,:), intent(inout) :: tcfc8
!--------------------------------------------------------------------
!--------------------------------------------------------------------
! intent(in)variables:
!
! Atmos_input
! tch4n2oe
!
! intent(inout) variables:
!
! Optical
! tcfc8
!
! intent(out) variables:
!
! emspec
! emspecf
!
!--------------------------------------------------------------------
!--------------------------------------------------------------------
! local variables
integer :: i,j,k,kk,m
real, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2), &
size(Atmos_input%temp,3)) :: &
temp, tflux, tpl1, tpl2, dte1, &
dte2, dxsp, ylog, dysp, emiss, &
emd1, emd2, dysp1
integer, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2), &
size(Atmos_input%temp,3)) :: &
ixsp, iysp, iysp1, ixoe1, ixoe2
real, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2), &
size(Atmos_input%temp,3)) :: ttmp
real, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2), &
size(Atmos_input%temp,3), NBTRGE) :: &
emissf, emd2f, emd1f
!--------------------------------------------------------------------
! local variables:
!
! temp
! tflux
! tpl1
! tpl2
! dte1
! dte2
! dxsp
! ylog
! dysp
! dysp1
! emiss
! emd1
! emd2
! ixsp
! iysp
! iysp1
! ixoe1
! ixoe2
! emissf
! emd2f
! emd1f
! m,k
!
!----------------------------------------------------------------------
!--------------------------------------------------------------------
!
!--------------------------------------------------------------------
do kk = 1,size(tflux(:,:,:),3)
do j = 1,size(tflux(:,:,:),2)
do i = 1,size(tflux(:,:,:),1)
tflux(i,j,kk) = Atmos_input%tflux(i,j,kk)
end do
end do
end do
do kk = 1,size(temp(:,:,:),3)
do j = 1,size(temp(:,:,:),2)
do i = 1,size(temp(:,:,:),1)
temp(i,j,kk) = Atmos_input%temp(i,j,kk)
end do
end do
end do
do j = 1,size(tpl1(:,:,:),2)
do i = 1,size(tpl1(:,:,:),1)
tpl1(i,j,KS) = temp(i,j,KE)
end do
end do
do kk = KS+1,KE
do j = 1,size(tpl1(:,:,:),2)
do i = 1,size(tpl1(:,:,:),1)
tpl1(i,j,kk) = tflux(i,j,kk)
end do
end do
end do
do j = 1,size(tpl1(:,:,:),2)
do i = 1,size(tpl1(:,:,:),1)
tpl1(i,j,KE+1) = 0.5E+00*(tflux(i,j,KE+1) + &
temp(i,j,KE))
end do
end do
do kk = KS+1,KE
do j = 1,size(tpl2(:,:,:),2)
do i = 1,size(tpl2(:,:,:),1)
tpl2(i,j,kk) = tflux(i,j,kk)
end do
end do
end do
do j = 1,size(tpl2(:,:,:),2)
do i = 1,size(tpl2(:,:,:),1)
tpl2(i,j,KE+1) = 0.5E+00*(tflux(i,j,KE) + &
temp(i,j,KE))
end do
end do
!--------------------------------------------------------------------
!
!--------------------------------------------------------------------
call locate_in_table(temp_1, temp, dte1, ixoe1, KS, KE+1)
call locate_in_table(temp_1, tflux, dte2, ixoe2, KS, KE+1)
ttmp(:,:,:)=ixoe2(:,:,:)
do kk = KS,KE
do j = 1,size(ixoe2(:,:,:),2)
do i = 1,size(ixoe2(:,:,:),1)
ixoe2(i,j,kk) = ttmp(i,j,kk+KS+1-(KS))
end do
end do
end do
ttmp(:,:,:)=dte2 (:,:,:)
do kk = KS,KE
do j = 1,size(dte2(:,:,:),2)
do i = 1,size(dte2(:,:,:),1)
dte2(i,j,kk) = ttmp(i,j,kk+KS+1-(KS))
end do
end do
end do
do j = 1,size(ixoe2(:,:,:),2)
do i = 1,size(ixoe2(:,:,:),1)
ixoe2(i,j,KE+1) = ixoe1(i,j,KE)
end do
end do
do j = 1,size(dte2(:,:,:),2)
do i = 1,size(dte2(:,:,:),1)
dte2(i,j,KE+1) = dte1 (i,j,KE)
end do
end do
do j = 1,size(ixsp(:,:,:),2)
do i = 1,size(ixsp(:,:,:),1)
ixsp(i,j,KE) = ixoe2(i,j,KE-1)
end do
end do
do j = 1,size(ixsp(:,:,:),2)
do i = 1,size(ixsp(:,:,:),1)
ixsp(i,j,KE+1) = ixoe1(i,j,KE-1)
end do
end do
do j = 1,size(dxsp(:,:,:),2)
do i = 1,size(dxsp(:,:,:),1)
dxsp(i,j,KE) = dte2(i,j,KE-1)
end do
end do
do j = 1,size(dxsp(:,:,:),2)
do i = 1,size(dxsp(:,:,:),1)
dxsp(i,j,KE+1) = dte1(i,j,KE-1)
end do
end do
if (Lw_control%do_h2o) then
do j = 1,size(ylog(:,:,:),2)
do i = 1,size(ylog(:,:,:),1)
ylog(i,j,KE ) = ALOG10(Optical%var2(i,j,KE))
end do
end do
do j = 1,size(ylog(:,:,:),2)
do i = 1,size(ylog(:,:,:),1)
ylog(i,j,KE+1) = ALOG10(Optical%var2(i,j,KE) + Optical%empl1(i,j,KE))
end do
end do
call locate_in_table (mass_1, ylog, dysp, iysp, KE, KE+1)
iysp(:,:,KE:KE+1) = iysp(:,:,KE:KE+1) + 1
!--------------------------------------------------------------------
! compute exchange terms in the flux equation for two terms used
! for nearby layer computations.
!--------------------------------------------------------------------
call looktab (tab2, ixsp, iysp, dxsp, dysp, emiss, KE, KE+1)
else
iysp1(:,:,KE:KE+1) = 1
dysp1(:,:,KE:KE+1) = 0.0
call looktab (tab2, ixsp, iysp1, dxsp, dysp1, emiss, KE, KE+1)
endif
!----------------------------------------------------------------------
! obtain index values of h2o pressure-scaled mass for each band
! in the 1200-1400 range.
!---------------------------------------------------------------------
if (nbtrge > 0) then
do m=1,NBTRGE
if (Lw_control%do_h2o) then
do j = 1,size(ylog(:,:,:),2)
do i = 1,size(ylog(:,:,:),1)
ylog(i,j,KE ) = ALOG10(Optical%vrpfh2o(i,j,KE,m))
end do
end do
do j = 1,size(ylog(:,:,:),2)
do i = 1,size(ylog(:,:,:),1)
ylog(i,j,KE+1) = ALOG10(Optical%vrpfh2o(i,j,KE,m) + Optical%empl1f(i,j,KE,m))
end do
end do
call locate_in_table (mass_1, ylog, dysp, iysp, KE, KE+1)
iysp(:,:,KE:KE+1) = iysp(:,:,KE:KE+1) + 1
!-----------------------------------------------------------------------
! compute exchange terms in the flux equation for two terms used
! for nearby layer computations.
!---------------------------------------------------------------------
call looktab (tab2a, ixsp, iysp, dxsp, dysp, &
emissf(:,:,:,m), KE, KE+1, m)
else
iysp1(:,:,KE:KE+1) = 1
dysp1(:,:,KE:KE+1) = 0.0
call looktab (tab2a, ixsp, iysp1, dxsp, dysp1, &
emissf(:,:,:,m), KE, KE+1, m)
endif
enddo
endif
!-----------------------------------------------------------------------
! compute nearby layer transmissivities for h2o.
!--------------------------------------------------------------------
if (Lw_control%do_h2o) then
call locate_in_table (temp_1, tpl1, dxsp, ixsp, KS, KE+1)
do kk = KS,KE+1
do j = 1,size(ylog(:,:,:),2)
do i = 1,size(ylog(:,:,:),1)
ylog(i,j,kk) = ALOG10(Optical%empl1(i,j,kk))
end do
end do
end do
call locate_in_table (mass_1, ylog, dysp, iysp, KS, KE+1)
iysp(:,:,KS:KE+1) = iysp(:,:,KS:KE+1) + 1
call looktab (tab3, ixsp, iysp, dxsp, dysp, emd1, KS, KE+1)
else
call locate_in_table (temp_1, tpl1, dxsp, ixsp, KS, KE+1)
iysp1(:,:,KS:KE+1) = 1
dysp1(:,:,KS:KE+1) = 0.0
call looktab (tab3, ixsp, iysp1, dxsp, dysp1, emd1, KS, KE+1)
endif
!----------------------------------------------------------------------
! obtain index values of h2o pressure-scaled mass for each band
! in the 1200-1400 range.
!------------------------------------------------------------------
if (nbtrge > 0) then
do m=1,NBTRGE
if (Lw_control%do_h2o) then
do kk = KS,KE+1
do j = 1,size(ylog(:,:,:),2)
do i = 1,size(ylog(:,:,:),1)
ylog(i,j,kk) = ALOG10(Optical%empl1f(i,j,kk,m))
end do
end do
end do
call locate_in_table (mass_1, ylog, dysp, iysp, KS, KE+1)
iysp(:,:,KS:KE+1) = iysp(:,:,KS:KE+1) + 1
call looktab (tab3a, ixsp, iysp, dxsp, dysp, &
emd1f(:,:,:,m), KS, KE+1, m)
else
iysp1(:,:,KS:KE+1) = 1
dysp1(:,:,KS:KE+1) = 0.0
call looktab (tab3a, ixsp, iysp1, dxsp, dysp1, &
emd1f(:,:,:,m), KS, KE+1, m)
endif
enddo
endif
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
if (Lw_control%do_h2o) then
call locate_in_table (temp_1, tpl2, dxsp, ixsp, KS+1, KE+1)
do kk = KS+1,KE+1
do j = 1,size(ylog(:,:,:),2)
do i = 1,size(ylog(:,:,:),1)
ylog(i,j,kk) = ALOG10(Optical%empl2(i,j,kk))
end do
end do
end do
call locate_in_table (mass_1, ylog, dysp, iysp, KS+1, KE+1)
iysp(:,:,KS+1:KE+1) = iysp(:,:,KS+1:KE+1) + 1
call looktab (tab3, ixsp, iysp, dxsp, dysp, emd2, KS+1, KE+1)
else
call locate_in_table (temp_1, tpl2, dxsp, ixsp, KS+1, KE+1)
iysp1(:,:,KS+1:KE+1) = 1
dysp1(:,:,KS+1:KE+1) = 0.0
call looktab (tab3, ixsp, iysp1, dxsp, dysp1, emd2, KS+1, KE+1)
endif
!----------------------------------------------------------------------
! obtain index values of h2o pressure-scaled mass for each band
! in the 1200-1400 range.
!---------------------------------------------------------------------
if (nbtrge > 0 ) then
do m=1,NBTRGE
if (Lw_control%do_h2o) then
do kk = KS+1,KE+1
do j = 1,size(ylog(:,:,:),2)
do i = 1,size(ylog(:,:,:),1)
ylog(i,j,kk) = ALOG10(Optical%empl2f(i,j,kk,m))
end do
end do
end do
call locate_in_table (mass_1, ylog, dysp, iysp, KS+1, KE+1)
iysp(:,:,KS+1:KE+1) = iysp(:,:,KS+1:KE+1) + 1
call looktab (tab3a, ixsp, iysp, dxsp, dysp, &
emd2f(:,:,:,m), KS+1, KE+1, m)
else
iysp1(:,:,KS+1:KE+1) = 1
dysp1(:,:,KS+1:KE+1) = 0.0
call looktab (tab3a, ixsp, iysp1, dxsp, dysp1, &
emd2f(:,:,:,m), KS+1, KE+1, m)
endif
enddo
endif
!----------------------------------------------------------------------
! compute nearby layer and special-case transmissivities for
! emissivity using methods for h2o given in reference (4).
!--------------------------------------------------------------------
if (Lw_control%do_h2o) then
do j = 1,size(emspec(:,:,:),2)
do i = 1,size(emspec(:,:,:),1)
emspec(i,j,KS ) = (emd1(i,j,KS)*Optical%empl1(i,j,KS) - &
emd1(i,j,KE+1)*Optical%empl1(i,j,KE+1))/ &
Optical%emx1(i,j) + 0.25E+00*(emiss(i,j,KE) + &
emiss(i,j,KE+1))
end do
end do
do j = 1,size(emspec(:,:,:),2)
do i = 1,size(emspec(:,:,:),1)
emspec(i,j,KS+1) = 2.0E+00*(emd1(i,j,KS)*Optical%empl1(i,j,KS) - &
emd2(i,j,KE+1)*Optical%empl2(i,j,KE+1))/ &
Optical%emx2(i,j)
end do
end do
else
do j = 1,size(emspec(:,:,:),2)
do i = 1,size(emspec(:,:,:),1)
emspec(i,j,KS ) = (emd1(i,j,KS) - emd1(i,j,KE+1)) / &
(Atmos_input%press(i,j,KE) - Atmos_input%pflux(i,j,KE)) + &
0.25E+00*(emiss(i,j,KE) + emiss(i,j,KE+1))
end do
end do
do j = 1,size(emspec(:,:,:),2)
do i = 1,size(emspec(:,:,:),1)
emspec(i,j,KS+1) = 2.0E+00*(emd1(i,j,KS) - emd2(i,j,KE+1)) / &
(Atmos_input%pflux(i,j,KE+1) - Atmos_input%press(i,j,KE))
end do
end do
endif
if (nbtrge > 0) then
do m=1,NBTRGE
if (Lw_control%do_h2o) then
do j = 1,size(emspecf(:,:,:,:),2)
do i = 1,size(emspecf(:,:,:,:),1)
emspecf(i,j,KS,m ) = (emd1f(i,j,KS,m)*Optical%empl1f(i,j,KS,m) - &
emd1f(i,j,KE+1,m)*Optical%empl1f(i,j,KE+1,m))/ &
Optical%emx1f(i,j,m) + 0.25E+00*(emissf(i,j,KE,m) + &
emissf(i,j,KE+1,m))
end do
end do
do j = 1,size(emspecf(:,:,:,:),2)
do i = 1,size(emspecf(:,:,:,:),1)
emspecf(i,j,KS+1,m) = 2.0E+00* &
(emd1f(i,j,KS,m)*Optical%empl1f(i,j,KS,m) - &
emd2f(i,j,KE+1,m)*Optical%empl2f(i,j,KE+1,m)) / &
Optical%emx2f(i,j,m)
end do
end do
else
do j = 1,size(emspecf(:,:,:,:),2)
do i = 1,size(emspecf(:,:,:,:),1)
emspecf(i,j,KS,m ) = (emd1f(i,j,KS,m) - emd1f(i,j,KE+1,m)) / &
(Atmos_input%press(i,j,KE) - Atmos_input%pflux(i,j,KE)) + &
0.25E+00*(emissf(i,j,KE,m) + emissf(i,j,KE+1,m))
end do
end do
do j = 1,size(emspecf(:,:,:,:),2)
do i = 1,size(emspecf(:,:,:,:),1)
emspecf(i,j,KS+1,m) = 2.0E+00* &
(emd1f(i,j,KS,m) - emd2f(i,j,KE+1,m)) / &
(Atmos_input%pflux(i,j,KE+1) - Atmos_input%press(i,j,KE))
end do
end do
endif
enddo
endif
!--------------------------------------------------------------------
! add the effects of other radiative gases on these flux arrays.
! the lbl transmissivities imply (at present) NBTRG = NBTRGE = 1).
! thus, tch4e and tn2oe are obtained directly from the transmission
! functions.
!----------------------------------------------------------------------
if (nbtrge > 0) then
if (Lw_control%do_ch4 .or. Lw_control%do_n2o) then
do k=KS,KS+1
do j = 1,size(emspecf(:,:,:,:),2)
do i = 1,size(emspecf(:,:,:,:),1)
emspecf(i,j,K,1) = emspecf(i,j,K,1)*tch4n2oe(i,j,KE+1,1)
end do
end do
end do
endif
!--------------------------------------------------------------------
! add cfc transmissivities if species which absorb in this fre-
! quency range are present.
!----------------------------------------------------------------------
if (Lw_control%do_cfc) then
call cfc_indx8_part (8, Optical, tcfc8, KE)
do k=KS,KS+1
do j = 1,size(emspecf(:,:,:,:),2)
do i = 1,size(emspecf(:,:,:,:),1)
emspecf(i,j,K,1) = emspecf(i,j,K,1)*tcfc8(i,j,KE+1)
end do
end do
end do
endif
endif ! (nbtrge > 0)
!------------------------------------------------------------------
end subroutine esfc
!######################################################################
!
!
! Subroutine to compute thermal layer emissivity using pre computed
! look up tables
!
!
! call enear (Atmos_input, emisdg, Optical, &
! emisdgf, tch4n2oe, tcfc8 )
!
!
! Atmospheric input data such as temperature and flux level temp
!
!
! Emissivity of thermal layers
!
!
! Optical depth of thermal layers
!
!
! Emissivity of thermal layers including effects of minor gas species
!
!
! CH4 and N2O transmission function
!
!
! CFC transmission function
!
!
!
subroutine enear (Atmos_input, emisdg, Optical, emisdgf, tch4n2oe, &
tcfc8)
!--------------------------------------------------------------------
!
!--------------------------------------------------------------------
type(atmos_input_type), intent(in) :: Atmos_input
real, dimension (:,:,:), intent(out) :: emisdg
type(optical_path_type), intent(inout) :: Optical
real, dimension (:,:,:,:), intent(out) :: emisdgf
real, dimension (:,:,:,:), intent(in) :: tch4n2oe
real, dimension (:,:,:), intent(inout) :: tcfc8
!--------------------------------------------------------------------
!--------------------------------------------------------------------
! intent(in) variables:
!
! Atmos_input
! tch4n2oe
!
! intent(inout) variables:
!
! Optical
! tcfc8
!
! intent(out) variables:
!
! emisdg
! emisdgf
!
!---------------------------------------------------------------------
!--------------------------------------------------------------------
! local variables:
integer :: i,j,kk,m
real, dimension (size(emisdg,1), &
size(emisdg,2), &
size(emisdg,3)) :: dte1, dte2
integer, dimension (size(emisdg,1), &
size(emisdg,2), &
size(emisdg,3)) :: ixoe1, ixoe2
real, dimension (size(emisdg,1), &
size(emisdg,2), &
size(emisdg,3)) :: ttmp
real, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2), &
size(Atmos_input%temp,3)) :: &
temp, tflux, tpl1, tpl2, &
dxsp, ylog, dysp, dysp1, emd1, emd2
integer, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2), &
size(Atmos_input%temp,3)) :: &
ixsp, iysp, iysp1
real, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2), &
size(Atmos_input%temp,3), NBTRGE) :: &
emd2f, emd1f
!--------------------------------------------------------------------
! local variables:
!
! dte1
! dte2
! ixoe1
! ixoe2
! temp
! tflux
! tpl1
! tpl2
! dxsp
! ylog
! dysp
! dysp1
! emiss
! emd1
! emd2
! ixsp
! iysp
! iysp1
! emissf
! emd2f
! emd1f
!
!--------------------------------------------------------------------
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
do kk = 1,size(tflux(:,:,:),3)
do j = 1,size(tflux(:,:,:),2)
do i = 1,size(tflux(:,:,:),1)
tflux(i,j,kk) = Atmos_input%tflux(i,j,kk)
end do
end do
end do
do kk = 1,size(temp(:,:,:),3)
do j = 1,size(temp(:,:,:),2)
do i = 1,size(temp(:,:,:),1)
temp(i,j,kk) = Atmos_input%temp(i,j,kk)
end do
end do
end do
do j = 1,size(tpl1(:,:,:),2)
do i = 1,size(tpl1(:,:,:),1)
tpl1(i,j,KS) = temp(i,j,KE)
end do
end do
do kk = KS+1,KE
do j = 1,size(tpl1(:,:,:),2)
do i = 1,size(tpl1(:,:,:),1)
tpl1(i,j,kk) = tflux(i,j,kk)
end do
end do
end do
do j = 1,size(tpl1(:,:,:),2)
do i = 1,size(tpl1(:,:,:),1)
tpl1(i,j,KE+1) = 0.5E+00*(tflux(i,j,KE+1) + &
temp(i,j,KE))
end do
end do
do kk = KS+1,KE
do j = 1,size(tpl2(:,:,:),2)
do i = 1,size(tpl2(:,:,:),1)
tpl2(i,j,kk) = tflux(i,j,kk)
end do
end do
end do
do j = 1,size(tpl2(:,:,:),2)
do i = 1,size(tpl2(:,:,:),1)
tpl2(i,j,KE+1) = 0.5E+00*(tflux(i,j,KE) + &
temp(i,j,KE))
end do
end do
!--------------------------------------------------------------------
!
!--------------------------------------------------------------------
call locate_in_table(temp_1, temp, dte1, ixoe1, KS, KE+1)
call locate_in_table(temp_1, tflux, dte2, ixoe2, KS, KE+1)
ttmp(:,:,:)=ixoe2(:,:,:)
do kk = KS,KE
do j = 1,size(ixoe2(:,:,:),2)
do i = 1,size(ixoe2(:,:,:),1)
ixoe2(i,j,kk) = ttmp(i,j,kk+KS+1-(KS))
end do
end do
end do
ttmp(:,:,:)=dte2 (:,:,:)
do kk = KS,KE
do j = 1,size(dte2(:,:,:),2)
do i = 1,size(dte2(:,:,:),1)
dte2(i,j,kk) = ttmp(i,j,kk+KS+1-(KS))
end do
end do
end do
do j = 1,size(ixoe2(:,:,:),2)
do i = 1,size(ixoe2(:,:,:),1)
ixoe2(i,j,KE+1) = ixoe1(i,j,KE)
end do
end do
do j = 1,size(dte2(:,:,:),2)
do i = 1,size(dte2(:,:,:),1)
dte2(i,j,KE+1) = dte1 (i,j,KE)
end do
end do
do j = 1,size(ixsp(:,:,:),2)
do i = 1,size(ixsp(:,:,:),1)
ixsp(i,j,KE) = ixoe2(i,j,KE-1)
end do
end do
do j = 1,size(ixsp(:,:,:),2)
do i = 1,size(ixsp(:,:,:),1)
ixsp(i,j,KE+1) = ixoe1(i,j,KE-1)
end do
end do
do j = 1,size(dxsp(:,:,:),2)
do i = 1,size(dxsp(:,:,:),1)
dxsp(i,j,KE) = dte2(i,j,KE-1)
end do
end do
do j = 1,size(dxsp(:,:,:),2)
do i = 1,size(dxsp(:,:,:),1)
dxsp(i,j,KE+1) = dte1(i,j,KE-1)
end do
end do
!--------------------------------------------------------------------
! compute exchange terms in the flux equation for two terms used
! for nearby layer computations.
! obtain index values of h2o pressure-scaled mass for each band
! in the 1200-1400 range.
!---------------------------------------------------------------------
!-----------------------------------------------------------------------
! compute nearby layer transmissivities for h2o.
!--------------------------------------------------------------------
if (Lw_control%do_h2o) then
call locate_in_table (temp_1, tpl1, dxsp, ixsp, KS, KE+1)
do kk = KS,KE+1
do j = 1,size(ylog(:,:,:),2)
do i = 1,size(ylog(:,:,:),1)
ylog(i,j,kk) = ALOG10(Optical%empl1(i,j,kk))
end do
end do
end do
call locate_in_table (mass_1, ylog, dysp, iysp, KS, KE+1)
iysp(:,:,KS:KE+1) = iysp(:,:,KS:KE+1) + 1
call looktab (tab3, ixsp, iysp, dxsp, dysp, emd1, KS, KE+1)
else
call locate_in_table (temp_1, tpl1, dxsp, ixsp, KS, KE+1)
iysp1(:,:,KS:KE+1) = 1
dysp1(:,:,KS:KE+1) = 0.0
call looktab (tab3, ixsp, iysp1, dxsp, dysp1, emd1, KS, KE+1)
endif
!----------------------------------------------------------------------
! obtain index values of h2o pressure-scaled mass for each band
! in the 1200-1400 range.
!------------------------------------------------------------------
if (nbtrge > 0) then
do m=1,NBTRGE
if (Lw_control%do_h2o) then
do kk = KS,KE+1
do j = 1,size(ylog(:,:,:),2)
do i = 1,size(ylog(:,:,:),1)
ylog(i,j,kk) = ALOG10(Optical%empl1f(i,j,kk,m))
end do
end do
end do
call locate_in_table (mass_1, ylog, dysp, iysp, KS, KE+1)
iysp(:,:,KS:KE+1) = iysp(:,:,KS:KE+1) + 1
call looktab (tab3a, ixsp, iysp, dxsp, dysp, &
emd1f(:,:,:,m), KS, KE+1, m)
else
iysp1(:,:,KS:KE+1) = 1
dysp1(:,:,KS:KE+1) = 0.0
call looktab (tab3a, ixsp, iysp1, dxsp, dysp1, &
emd1f(:,:,:,m), KS, KE+1, m)
endif
enddo
endif
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
if (Lw_control%do_h2o) then
call locate_in_table (temp_1, tpl2, dxsp, ixsp, KS+1, KE+1)
do kk = KS+1,KE+1
do j = 1,size(ylog(:,:,:),2)
do i = 1,size(ylog(:,:,:),1)
ylog(i,j,kk) = ALOG10(Optical%empl2(i,j,kk))
end do
end do
end do
call locate_in_table (mass_1, ylog, dysp, iysp, KS+1, KE+1)
iysp(:,:,KS+1:KE+1) = iysp(:,:,KS+1:KE+1) + 1
call looktab (tab3, ixsp, iysp, dxsp, dysp, emd2, KS+1, KE+1)
else
call locate_in_table (temp_1, tpl2, dxsp, ixsp, KS+1, KE+1)
iysp1(:,:,KS+1:KE+1) = 1
dysp1(:,:,KS+1:KE+1) = 0.0
call looktab (tab3, ixsp, iysp1, dxsp, dysp1, emd2, KS+1, KE+1)
endif
!----------------------------------------------------------------------
! obtain index values of h2o pressure-scaled mass for each band
! in the 1200-1400 range.
!---------------------------------------------------------------------
if (nbtrge > 0) then
do m=1,NBTRGE
if (Lw_control%do_h2o) then
do kk = KS+1,KE+1
do j = 1,size(ylog(:,:,:),2)
do i = 1,size(ylog(:,:,:),1)
ylog(i,j,kk) = ALOG10(Optical%empl2f(i,j,kk,m))
end do
end do
end do
call locate_in_table (mass_1, ylog, dysp, iysp, KS+1, KE+1)
iysp(:,:,KS+1:KE+1) = iysp(:,:,KS+1:KE+1) + 1
call looktab (tab3a, ixsp, iysp, dxsp, dysp, &
emd2f(:,:,:,m), KS+1, KE+1, m)
else
iysp1(:,:,KS+1:KE+1) = 1
dysp1(:,:,KS+1:KE+1) = 0.0
call looktab (tab3a, ixsp, iysp1, dxsp, dysp1, &
emd2f(:,:,:,m), KS+1, KE+1, m)
endif
enddo
endif
!----------------------------------------------------------------------
! compute nearby layer and special-case transmissivities for
! emissivity using methods for h2o given in reference (4).
!--------------------------------------------------------------------
do kk = KS+1,KE
do j = 1,size(emisdg(:,:,:),2)
do i = 1,size(emisdg(:,:,:),1)
emisdg(i,j,kk) = emd2(i,j,kk) + emd1(i,j,kk)
end do
end do
end do
do j = 1,size(emisdg(:,:,:),2)
do i = 1,size(emisdg(:,:,:),1)
emisdg(i,j,KE+1) = 2.0E+00*emd1(i,j,KE+1)
end do
end do
if (nbtrge > 0) then
do m=1,NBTRGE
do kk = KS+1,KE
do j = 1,size(emisdgf(:,:,:,:),2)
do i = 1,size(emisdgf(:,:,:,:),1)
emisdgf(i,j,kk,m) = &
emd2f(i,j,kk,m) + emd1f(i,j,kk,m)
end do
end do
end do
do j = 1,size(emisdgf(:,:,:,:),2)
do i = 1,size(emisdgf(:,:,:,:),1)
emisdgf(i,j,KE+1,m) = 2.0E+00*emd1f(i,j,KE+1,m)
end do
end do
enddo
endif
!--------------------------------------------------------------------
! add the effects of other radiative gases on these flux arrays.
! the lbl transmissivities imply (at present) NBTRG = NBTRGE = 1).
! thus, tch4e and tn2oe are obtained directly from the transmission
! functions.
!----------------------------------------------------------------------
if (nbtrge > 0) then
if (Lw_control%do_ch4 .or. Lw_control%do_n2o) then
do kk = KS+1,KE+1
do j = 1,size(emisdgf(:,:,:,:),2)
do i = 1,size(emisdgf(:,:,:,:),1)
emisdgf(i,j,kk,1) = emisdgf(i,j,kk,1) * &
tch4n2oe(i,j,kk,1)
end do
end do
end do
endif
!--------------------------------------------------------------------
! add cfc transmissivities if species which absorb in this fre-
! quency range are present.
!----------------------------------------------------------------------
if (Lw_control%do_cfc) then
call cfc_indx8_part (8, Optical, tcfc8, KE)
do kk = KS+1,KE+1
do j = 1,size(emisdgf(:,:,:,:),2)
do i = 1,size(emisdgf(:,:,:,:),1)
emisdgf(i,j,kk,1) = emisdgf(i,j,kk,1) * &
tcfc8(i,j,kk)
end do
end do
end do
endif
endif ! (nbtrge > 0)
!------------------------------------------------------------------
end subroutine enear
!####################################################################
!
!
! Subroutine to calculate CO2 source function
!
!
! call co2_source_calc (Atmos_input, Rad_gases, sorc, Gas_tf, &
! source_band, dsrcdp_band)
!
!
! Atmospheric input data
!
!
! Radiative gases properties
!
!
! CO2 source function results
!
!
! Gas transmission functions
!
!
! CO2 source function bands
!
!
! Difference of source function between nearby thermal layers
!
!
!
subroutine co2_source_calc (Atmos_input, Rad_gases, sorc, Gas_tf, &
source_band, dsrcdp_band)
!--------------------------------------------------------------------
!
!--------------------------------------------------------------------
type(atmos_input_type), intent(in) :: Atmos_input
type(radiative_gases_type), intent(in) :: Rad_gases
real, dimension(:,:,:,:), intent(out) :: sorc
type(gas_tf_type), intent(in) :: Gas_tf
real, dimension(:,:,:,:), intent(out) :: source_band, dsrcdp_band
!----------------------------------------------------------------------
! intent(in) variables:
!
! Atmos_input
! Rad_gases
! Gas_tf
!
! intent(out) variables:
!
! sorc
! source_band
! dsrcdp_band
!
!----------------------------------------------------------------------
!---------------------------------------------------------------------
! local variables:
real, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2), &
size(Atmos_input%temp,3)) :: &
dte1, press, pflux, temp
integer, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2), &
size(Atmos_input%temp,3)) :: ixoe1
integer :: i,j,kk
integer :: n, ioffset, m
integer :: nbly
real :: rrvco2
!---------------------------------------------------------------------
! local variables:
!
! dte1
! press
! pflux
! temp
! ixoe1
! n
! ioffset
! m
! nbly
! rrvco2
!
!---------------------------------------------------------------------
!--------------------------------------------------------------------
!
!--------------------------------------------------------------------
ioffset = Lw_parameters%offset
nbly = 16+ioffset
!--------------------------------------------------------------------
! convert press and pflux to cgs.
do kk = 1,size(press(:,:,:),3)
do j = 1,size(press(:,:,:),2)
do i = 1,size(press(:,:,:),1)
press(i,j,kk) = 10.0*Atmos_input%press(i,j,kk)
end do
end do
end do
do kk = 1,size(pflux(:,:,:),3)
do j = 1,size(pflux(:,:,:),2)
do i = 1,size(pflux(:,:,:),1)
pflux(i,j,kk) = 10.0*Atmos_input%pflux(i,j,kk)
end do
end do
end do
do kk = 1,size(temp(:,:,:),3)
do j = 1,size(temp(:,:,:),2)
do i = 1,size(temp(:,:,:),1)
temp(i,j,kk) = Atmos_input%temp(i,j,kk)
end do
end do
end do
rrvco2 = Rad_gases%rrvco2
!----------------------------------------------------------------------
! compute source function for frequency bands (9+ioffset to NBLY-1)
! at layer temperatures using table lookup.
!----------------------------------------------------------------------
call locate_in_table(temp_1, temp, dte1, ixoe1, KS, Ke+1)
do n=9+ioffset,NBLY-1
call looktab (tabsr, ixoe1, dte1, &
sorc(:,:,:,n-ioffset-8), KS, ke+1, n)
enddo
!-----------------------------------------------------------------------
! compute the nlte source function for co2.
!-----------------------------------------------------------------------
if (do_nlte) then
call nlte (pflux, press, rrvco2, sorc, Gas_tf)
endif
!----------------------------------------------------------------------
! define "source function" appropriate for emissivity calculations
! (temp**4), source functions for selected ranges including more
! than 1 frequency band (sorc15 for 15 um co2 band)
! and differences in source functions (deltab) over
! pressure layers.
!
! note: the values of sorc, sorc15, sorcwin, and derivatives
! depend on the no. of freq. bands!
!-----------------------------------------------------------------------
do kk = 1,size(source_band(:,:,:,:),3)
do j = 1,size(source_band(:,:,:,:),2)
do i = 1,size(source_band(:,:,:,:),1)
source_band(i,j,kk,1) = Atmos_input%temp (i,j,kk+KS-(1))**4
end do
end do
end do
do kk = 1,size(source_band(:,:,:,:),3)
do j = 1,size(source_band(:,:,:,:),2)
do i = 1,size(source_band(:,:,:,:),1)
source_band(i,j,kk,2) = sorc(i,j,kk, 1) + &
sorc(i,j,kk, 2) + &
sorc(i,j,kk, 3 )
end do
end do
end do
do kk = 1,size(source_band(:,:,:,:),3)
do j = 1,size(source_band(:,:,:,:),2)
do i = 1,size(source_band(:,:,:,:),1)
source_band(i,j,kk,3) = sorc(i,j,kk,4 )
end do
end do
end do
do kk = 1,size(source_band(:,:,:,:),3)
do j = 1,size(source_band(:,:,:,:),2)
do i = 1,size(source_band(:,:,:,:),1)
source_band(i,j,kk,4) = sorc(i,j,kk,5 )
end do
end do
end do
do kk = 1,size(source_band(:,:,:,:),3)
do j = 1,size(source_band(:,:,:,:),2)
do i = 1,size(source_band(:,:,:,:),1)
source_band(i,j,kk,5) = sorc(i,j,kk, 6 )
end do
end do
end do
do kk = 1,size(source_band(:,:,:,:),3)
do j = 1,size(source_band(:,:,:,:),2)
do i = 1,size(source_band(:,:,:,:),1)
source_band(i,j,kk,6) = sorc(i,j,kk,7 )
end do
end do
end do
do m=1,NBTRGE
do kk = 1,size(source_band(:,:,:,:),3)
do j = 1,size(source_band(:,:,:,:),2)
do i = 1,size(source_band(:,:,:,:),1)
source_band(i,j,kk,6+m) = source_band(i,j,kk,1)
end do
end do
end do
end do
do n=1, 6+NBTRGE
do kk = KS+1,KE+1
do j = 1,size(dsrcdp_band(:,:,:,:),2)
do i = 1,size(dsrcdp_band(:,:,:,:),1)
dsrcdp_band(i,j,kk,n) = source_band(i,j,kk,n) - &
source_band(i,j,kk+KS-(KS+1),n)
end do
end do
end do
end do
!-------------------------------------------------------------------
end subroutine co2_source_calc
!#####################################################################
!
!
! nlte is the present formulation of an nlte calculation of the
! source function in the 15 um region (two bands).
!
!
! nlte is the present formulation of an nlte calculation of the
! source function in the 15 um region (two bands).
!
! the essential theory is:
!
! phi = C*j
! j = b + E*phi
!
! where
! C = Curtis matrix
! E = NLTE contribution (diagonal matrix)
! phi = heating rate vector
! b = LTE source function vector
! j = NLTE source function vector
!
! j = b (by assumption) for pressure layers > ixnltr
! j = b (by assumption) for pressure layers > ixprnlte
! E is obtained using a formulation devised by Fels (denoted
! Ri in his notes).
!
!
! call nlte (pflux, press, rrvco2, sorc, Gas_tf)
!
!
! pressure values at flux levels.
!
!
! pressure cordinates
!
!
! CO2 volumn mixing ratio
!
!
! CO2 source function to be calculated
!
!
! Gas transmission function
!
!
!
subroutine nlte (pflux, press, rrvco2, sorc, Gas_tf)
!-----------------------------------------------------------------------
! nlte is the present formulation of an nlte calculation of the
! source function in the 15 um region (two bands).
!
! the essential theory is:
!
! phi = C*j
! j = b + E*phi
!
! where
! C = Curtis matrix
! E = NLTE contribution (diagonal matrix)
! phi = heating rate vector
! b = LTE source function vector
! j = NLTE source function vector
!
! j = b (by assumption) for pressure layers > ixnltr
! j = b (by assumption) for pressure layers > ixprnlte
! E is obtained using a formulation devised by Fels (denoted
! Ri in his notes).
!
! author: m. d. schwarzkopf
!
! revised: 1/1/93
!
! certified: radiation version 1.0
!-----------------------------------------------------------------------
real, dimension (:,:,:), intent(in) :: pflux, press
real, intent(in) :: rrvco2
real, dimension (:,:,:,:), intent(inout) :: sorc
type(gas_tf_type), intent(in) :: Gas_tf
!---------------------------------------------------------------------
! intent(in) variables:
!
! pflux
! press
! rrvco2
! Gas_tf
!
! intent(inout) variables:
!
! sorc
!
!--------------------------------------------------------------------
!---------------------------------------------------------------------
! local variables:
real, dimension (size(press,1), size(press,2), &
ixprnlte) :: &
ag, az, bdenom, cdiag, &
tcoll, phifx, phivar
real, dimension (size(press,1), size(press,2), &
ixprnlte, NBCO215) :: &
fnlte
real, dimension (size(press,1), size(press,2), &
size(press,3)-1, ixprnlte ) :: &
cmtrx
real :: degen = 0.5
integer :: i,j,kk
integer :: n, k, inb, kp, ioffset
!---------------------------------------------------------------------
! local variables:
!
! ag
! az
! bdenom
! cdiag
! tcoll
! phifx fixed portion of PHI (contributions from
! layers > ixnltr, where j(k) = b(k))
! layers > ixprnlte, where j(k) = b(k))
! phivar varying portion of PHI (contributions
! from layers <= ixprnlte).
! from layers <= ixnltr).
! fnlte NLTE contribution: (E in above notes)
! cmtrx
! degen degeneracy factor (= 0.5)
! n
! k
! inb
! kp
! ioffset
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!---------------------------------------------------------------------
ioffset = Lw_parameters%offset
!--------------------------------------------------------------------
!-----------------------------------------------------------------------
! compute curtis matrix for both frequency bands.
!-----------------------------------------------------------------------
call co2curt (pflux, cmtrx, Gas_tf)
do k=KS,ixprnlte
do j = 1,size(cdiag(:,:,:),2)
do i = 1,size(cdiag(:,:,:),1)
cdiag(i,j,k) = cmtrx(i,j,k,k)
end do
end do
end do
!-----------------------------------------------------------------------
! collisional relaxation time (see fels notes for "tcoll")
!-----------------------------------------------------------------------
do k=KS,ixprnlte
do j = 1,size(tcoll(:,:,:),2)
do i = 1,size(tcoll(:,:,:),1)
tcoll(i,j,k) = degen*1.5E-05*press(i,j,KE+1)/ &
(seconds_per_day*press(i,j,k))
end do
end do
end do
!-----------------------------------------------------------------------
! compute NLTE contribution for each band at each pressure level
! <= ixprnlte. fnlte = zero by assumption at other levels.
!-----------------------------------------------------------------------
do n=1,NBCO215
do kk = KS,ixprnlte
do j = 1,size(fnlte(:,:,:,:),2)
do i = 1,size(fnlte(:,:,:,:),1)
fnlte(i,j,kk,n) = 3.5E+00*tcoll(i,j,kk)* &
c1b7(n)/(rrvco2*c2b7(n))
end do
end do
end do
enddo
!-----------------------------------------------------------------------
! begin computations for (NBCO215) bands in 15um range.
!-----------------------------------------------------------------------
do inb = 1,NBCO215
do kk = KS,ixprnlte
do j = 1,size(bdenom(:,:,:),2)
do i = 1,size(bdenom(:,:,:),1)
bdenom(i,j,kk) = 1.0E+00/ &
(1.0E+00 - fnlte(i,j,kk,inb)* &
cdiag(i,j,kk))
end do
end do
end do
do kk = KS,ixprnlte
do j = 1,size(phifx(:,:,:),2)
do i = 1,size(phifx(:,:,:),1)
phifx(i,j,kk) = 0.0E+00
end do
end do
end do
do k=KS,ixprnlte
do kp=ixprnlte+1,KE
do j = 1,size(phifx(:,:,:),2)
do i = 1,size(phifx(:,:,:),1)
phifx(i,j,k) = phifx(i,j,k) + &
cmtrx(i,j,kp,k)*sorc(i,j,kp,inb )
end do
end do
end do
end do
do kk = KS,ixprnlte
do j = 1,size(az(:,:,:),2)
do i = 1,size(az(:,:,:),1)
az(i,j,kk) = sorc (i,j,kk,inb ) + &
fnlte(i,j,kk,inb)*phifx(i,j,kk)
end do
end do
end do
!----------------------------------------------------------------------
! first iteration. (J(k) = B(k)) as initial guess)
!-----------------------------------------------------------------------
do kk = KS,ixprnlte
do j = 1,size(phivar(:,:,:),2)
do i = 1,size(phivar(:,:,:),1)
phivar(i,j,kk) = 0.0E+00
end do
end do
end do
do k=KS,ixprnlte
do kp=KS,ixprnlte
do j = 1,size(phivar(:,:,:),2)
do i = 1,size(phivar(:,:,:),1)
phivar(i,j,k) = phivar(i,j,k) + &
cmtrx(i,j,kp,k)*sorc(i,j,kp,inb )
end do
end do
end do
end do
do kk = KS,ixprnlte
do j = 1,size(ag(:,:,:),2)
do i = 1,size(ag(:,:,:),1)
ag(i,j,kk) = fnlte(i,j,kk,inb)* &
(phivar(i,j,kk) - &
cdiag(i,j,kk)* &
sorc(i,j,kk,inb ))
end do
end do
end do
do kk = KS,ixprnlte
do j = 1,size(sorc(:,:,:,:),2)
do i = 1,size(sorc(:,:,:,:),1)
sorc(i,j,kk,inb ) = bdenom(i,j,kk)*&
(az(i,j,kk) + &
ag(i,j,kk))
end do
end do
end do
!-----------------------------------------------------------------------
! second iteration. (J(k) = result of first iteration as guess)
!-----------------------------------------------------------------------
do kk = KS,ixprnlte
do j = 1,size(phivar(:,:,:),2)
do i = 1,size(phivar(:,:,:),1)
phivar(i,j,kk) = 0.0E+00
end do
end do
end do
do k=KS,ixprnlte
do kp=KS,ixprnlte
do j = 1,size(phivar(:,:,:),2)
do i = 1,size(phivar(:,:,:),1)
phivar(i,j,k) = phivar(i,j,k) + &
cmtrx(i,j,kp,k)*sorc(i,j,kp,inb )
end do
end do
end do
end do
do kk = KS,ixprnlte
do j = 1,size(ag(:,:,:),2)
do i = 1,size(ag(:,:,:),1)
ag(i,j,kk) = fnlte(i,j,kk,inb)* &
(phivar(i,j,kk) - &
cdiag(i,j,kk)*sorc(i,j,kk,inb ))
end do
end do
end do
do kk = KS,ixprnlte
do j = 1,size(sorc(:,:,:,:),2)
do i = 1,size(sorc(:,:,:,:),1)
sorc(i,j,kk,inb ) = bdenom(i,j,kk)*&
(az(i,j,kk) + &
ag(i,j,kk))
end do
end do
end do
enddo
!-----------------------------------------------------------------------
end subroutine nlte
!#####################################################################
!
!
! co2curt computes Curtis matrix elements derived from co2
! transmission functions.
!
!
! call co2curt (pflux, cmtrx, Gas_tf)
!
!
! pressure values at flux levels
!
!
! Curtis matrix elements
!
!
! gas transmission function
!
!
!
subroutine co2curt (pflux, cmtrx, Gas_tf)
!----------------------------------------------------------------------
! co2curt computes Curtis matrix elements derived from co2
! transmission functions.
! functions.
!
! author: m. d. schwarzkopf
!
! revised: 8/18/94
!
! certified: radiation version 1.0
!
!---------------------------------------------------------------------
real, dimension(:,:, :), intent(in) :: pflux
real, dimension(:,:, :,:), intent(out) :: cmtrx
type(gas_tf_type), intent(in) :: Gas_tf
!---------------------------------------------------------------------
! intent(in) variables:
!
! pflux
! Gas_tf
!
! intent(out) variables:
!
! cmtrx cutris matrix.
!
!---------------------------------------------------------------------
!---------------------------------------------------------------------
! local variables:
real, dimension (size(pflux,1), &
size(pflux,2), &
size(pflux,3)-1) :: pdfinv
real, dimension (size(pflux,1), &
size(pflux,2), &
size(pflux,3)) :: co2row, co2rowp
integer :: k, krow, kp
integer :: i, j, kk
!---------------------------------------------------------------------
! local variables:
!
! pdfinv
! co2row
! co2rowp
! k
! krow
! kp
!
!---------------------------------------------------------------------
!-----------------------------------------------------------------------
! compute co2 transmission functions.
!-----------------------------------------------------------------------
do kk = KS,KE+1
do j = 1,size(co2row(:,:,:),2)
do i = 1,size(co2row(:,:,:),1)
co2row(i,j,kk) = 1.0E+00
end do
end do
end do
do kk = KS,KE+1
do j = 1,size(co2rowp(:,:,:),2)
do i = 1,size(co2rowp(:,:,:),1)
co2rowp(i,j,kk) = 1.0E+00
end do
end do
end do
!-----------------------------------------------------------------------
! compute curtis matrix for rows from KS to ixprnlte
!-----------------------------------------------------------------------
do k = KS,ixprnlte
krow = k
do j = 1,size(pdfinv(:,:,:),2)
do i = 1,size(pdfinv(:,:,:),1)
pdfinv(i,j,k) = 1.0/(pflux(i,j,k+1) - pflux(i,j,k))
end do
end do
call transcol ( KS, krow, KS, KE+1, co2row, Gas_tf)
call transcol ( KS, krow+1, KS, KE+1, co2rowp, Gas_tf)
do kp=KS,KE-1
do j = 1,size(cmtrx(:,:,:,:),2)
do i = 1,size(cmtrx(:,:,:,:),1)
cmtrx(i,j,kp,k) = radcon*pdfinv(i,j,k)* &
(co2rowp(i,j,kp) - co2rowp(i,j,kp+1) - &
co2row(i,j,kp) + co2row(i,j,kp+1))
end do
end do
end do
do j = 1,size(cmtrx(:,:,:,:),2)
do i = 1,size(cmtrx(:,:,:,:),1)
cmtrx(i,j,KE,k) = radcon*pdfinv(i,j,k)* &
(co2rowp(i,j,KE) - co2row(i,j,KE))
end do
end do
enddo
!--------------------------------------------------------------------
end subroutine co2curt
!####################################################################
!
!
! Calculate CO2 absorption coefficient based on its volume
! mixing ratio using precomputed lbl tables
!
!
! Calculate CO2 absorption coefficient based on its volume
! mixing ratio using precomputed lbl tables
!
!
! call co2_time_vary ( rrvco2 )
!
!
! CO2 volume mixing ratio
!
!
!
subroutine co2_time_vary ( rrvco2 )
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
real, intent(in ) :: rrvco2
!---------------------------------------------------------------------
! intent(in) variables:
!
! rrvco2
!
!---------------------------------------------------------------------
!----------------------------------------------------------------------
! local variables:
real :: co2_vmr !
!--------------------------------------------------------------------
!
!--------------------------------------------------------------------
co2_vmr = rrvco2*1.0E+06
!--------------------------------------------------------------------
!
!--------------------------------------------------------------------
call co2_lblinterp (co2_vmr )
!--------------------------------------------------------------------
end subroutine co2_time_vary
!####################################################################
!
!
! Calculate CH4 and N2O absorption coefficients from their
! mixing ratios using precomputed lbl tables
!
!
! Calculate CH4 and N2O absorption coefficients from their
! mixing ratios using precomputed lbl tables
!
!
! call ch4_n2o_time_vary (rrvch4, rrvn2o)
!
!
! ch4 volume mixing ratio
!
!
! n2o volume mixing ratio
!
!
!
subroutine ch4_time_vary (rrvch4)
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
real, intent(in) :: rrvch4
!---------------------------------------------------------------------
! intent(in) variables:
!
! rrvch4
!
!---------------------------------------------------------------------
!---------------------------------------------------------------------
! local variables:
real :: ch4_vmr !
!---------------------------------------------------------------------
! the ch4 volume mixing ratio is set to the initial value (rch4) and
! the mass mixing ratio is defined on the first access of this
! routine. then the lbl transmission function is calculated. after
! first access, this routine does nothing.
!--------------------------------------------------------------------
ch4_vmr = rrvch4*1.0E+09
call Ch4_lblinterp (ch4_vmr)
!---------------------------------------------------------------------
! the n2o volume mixing ratio is set to initial value (rn2o) and the
! mass mixing ratio is defined on the first access of this routine.
! routines are called to calculate the lbl transmission functions for
! n2o. after first access, this routine does nothing.
!--------------------------------------------------------------------
! n2o_vmr = rrvn2o*1.0E+09
! call N2o_lblinterp (n2o_vmr)
!----------------------------------------------------------------------
end subroutine ch4_time_vary
!####################################################################
!
!
! Calculate CH4 and N2O absorption coefficients from their
! mixing ratios using precomputed lbl tables
!
!
! Calculate CH4 and N2O absorption coefficients from their
! mixing ratios using precomputed lbl tables
!
!
! call ch4_n2o_time_vary (rrvch4, rrvn2o)
!
!
! ch4 volume mixing ratio
!
!
! n2o volume mixing ratio
!
!
!
subroutine n2o_time_vary (rrvn2o)
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
real, intent(in) :: rrvn2o
!---------------------------------------------------------------------
! intent(in) variables:
!
! rrvch4
! rrvn2o
!
!---------------------------------------------------------------------
!---------------------------------------------------------------------
! local variables:
! real :: ch4_vmr !
real :: n2o_vmr !
!---------------------------------------------------------------------
! the ch4 volume mixing ratio is set to the initial value (rch4) and
! the mass mixing ratio is defined on the first access of this
! routine. then the lbl transmission function is calculated. after
! first access, this routine does nothing.
!--------------------------------------------------------------------
! ch4_vmr = rrvch4*1.0E+09
! call Ch4_lblinterp (ch4_vmr)
!---------------------------------------------------------------------
! the n2o volume mixing ratio is set to initial value (rn2o) and the
! mass mixing ratio is defined on the first access of this routine.
! routines are called to calculate the lbl transmission functions for
! n2o. after first access, this routine does nothing.
!--------------------------------------------------------------------
n2o_vmr = rrvn2o*1.0E+09
call N2o_lblinterp (n2o_vmr)
!----------------------------------------------------------------------
end subroutine n2o_time_vary
!####################################################################
end module sealw99_mod